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ABSTRACT 


Two  approximate  parametric  interval  estimation  methods  for  system  reliability  using 
component  test  data  are  developed  and  evaluated.  One  method  can  be  applied  to  any 
coherent  system  with  components  which  have  exponential  failure  times  with  possibly 
different  failure  rates  and  different  mission  operating  times.  This  method  estimates  the 
ratios  of  component  failure  rates  which  are  then  used  to  develop  the  approximate  lower 
confidence  limit.  These  ratio  estimates  are  developed  with  and  without  jacknife  methods 
and  the  two  results  are  compared.  This  procedure  is  very  accurate  and  simple  to  com¬ 
pute,  requiring  the  use  of  standard  chi-square  tables.  This  ratio  method  is  subsequently 
extended  to  coherent  systems  with  components  whose  failure  times  have  a  Weibull  dis¬ 
tribution.  A  nearly  exact  parametric  lower  confidence  limit  for  P(X  >  *)  is  developed  and 
evaluated  where  .t  is  given  and  X  has  a  normal  distribution  with  unknown  mean  and 
variance.  This  procedure  is  also  simple  to  evaluate  and  requires  the  use  of  Student  t  ta¬ 
bles. 
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I.  BACKGROUND 


Parametric  confidence  interval  procedures  for  the  reliability  of  mechanical  systems 
are  much  less  developed  than  procedures  for  electrical  systems.  This  is  due  to  the  more 
complicated  failure  distributions  used  to  model  mechanical  hardware.  The  failure  rate 
of  a  senes  system  of  independent  components,  with  exponentially  distributed  failure 
times  and  equal  mission  times,  is  the  sum  of  the  failure  component  failure  rates.  This 
proper';;  has  permitted  development  of  numerous  methods  for  system  reliability  of  series 
and  other  coherent  systems  using  component  failure  data.  The  Weibull  and  extreme 
value  distributions  have  been  used  in  life  testing  methods  for  both  electrical  and  me¬ 
chanical  devices.  Several  methods  have  been  used  for  obtaining  point  estimates  of  pa¬ 
rameters  for  these  distributions  and  for  the  reliability  function  itself,  Harter  and  Moore 
[Ref.  1:  Pp.  889-901],  Mann  [Ref.  2:  pp.  231-256]. 

The  derivation  of  simple  confidence  limits  for  the  reliability  function  for  the  extreme 
value  distribution  with  parameters  z  and  <5  has  posed  probiems,  because  methods  based 

A 

on  a  pivotal  quantity  such  as  (i  -  z)IS  are  inadequate.  Lawless  [Ref.  3:  pp.  355-364], 
Johns  and  Lieberman  [Ref.  4:  pp.  135-175]  and  Thoman,  Bain  and  Antle  [Ref.  5:  pp. 
363-372]  have  developed  nearly  exact  procedures  for  confidence  limits  for  the  reliability 
function  of  the  Weibull  and  extreme  value  distributions.  Schneider  and  Weissfeld  [Ref. 
6:  pp.  179-186]  have  developed  interval  estimation  methods  for  percentiles  of  the  Weibull 
and  extreme  value  distribution  based  on  censored  data.  Although  extensive  methods 
have  been  developed  for  interval  estimates  of  the  reliability  of  a  single  component  with 
the  Weibull  or  extreme  value  failure  distribution,  very  few  parametric  interval  methods 
have  been  developed  for  system  reliability  using  component  test  data  with  Weibull  failure 
distributions. 

Two  approximate  parametric  interval  estimation  methods  for  reliability  of  coherent 
systems  using  component  test  data  are  developed  and  evaluated  in  this  thesis.  Evalu¬ 
ations  of  these  procedures  are  performed  using  computer  simulation  for  series  systems 
only.  One  of  these  methods  is  developed  for  the  reliability  of  a  series  system  whose 
components  have  exponential  failure  distributions  and  different  mission  times.  This 
procedure  was  found  to  be  quite  accurate  and  can  be  applied  to  coherent  systems  in 
general.  This  first  procedure  is  then  extended  to  the  case  where  components  of  the  sys- 


1 


tem  have  Weibull  failure  distributions.  The  method  used  is  an  extension  of  a  non- 
paramctric  method  developed  by  Myhri,  J.,  Rosenfeld,  A.  and  Saunders,  S.  [Ref.  7: 
pp.213-227]. 

The  normal  distribution  is  used  extensively  in  some  mechanical  reliabilty  models. 
Maximum  likelihood  and  minimum  variance  unbiased  estimators  for  P{X^.x0)  when 
both  n  and  o 2  are  unknown  were  developed  over  thirty  years  ago,  Lieberman  and 
Resnikoff  [Ref.  8:  pp.  457-516],  Folks  and  others  [Ref.  9:  pp.  43-50],  and  Barton  [Ref. 
10:  pp.  227-229].  Exact  interval  estimation  procedures  for  P{X>x0)  were  developed  by 
Owen  and  Hua  [Ref.  11:  pp.285-311]  using  the  non-central ;  distribution.  Letting  X  de¬ 
note  strength  and  Y  denote  stress,  mechanical  reliability  is  sometimes  modeled  as  the 
value  for  P(X  >  F).  Approximate  interval  estimation  procedures  for  this  probability 
when  X  and  Y  are  assumed  to  be  normal  have  been  developed  by  Church  and  Harris 
[Ref.  12:  pp.  49-54]  when  the  mean  and  variance  of  Fare  known.  Downton  [Ref.  13:  pp. 
551-55S]  modifies  their  procedure  slightly  to  get  more  accurate  bounds  and  suggests  an 
approximate  procedure  when  the  means  and  variances  of  both  X  and  Y  are  unknown. 
Lower  confidence  intervals  for  P{X  >  Y)  obtained  under  the  assumption  of  normality  for 
X  and  Y  can  lead  to  serious  error  when  either  X  or  Y  or  both  are  truncated  well  into  the 
tails.  Consequently,  P[X  >  xe)  may  be  a  more  reasonable  model  of  mechanical  reliability 
where  .v0  is  chosen  conservatively.  A  very  accurate  approximate  lower  confidence  limit 
procedure  for  P(X^x0)  is  developed  and  evaluated  in  this  thesis.  It  can  be  computed 
easily. 


II.  INTERVAL  ESTIMATION  PROCEDURE  -  EXPONENTIAL  CASE 


A  system  of  independent  components  is  coherent  if  an  increase  in  reliability  of  any 
one  of  its  components  does  not  cause  a  degredation  in  system  reliability.  Suppose  a 
coherent  system  has  k  components.  We  assume  that  the  failure  distribution  of  compo¬ 
nent  i  is  exponential  with  failure  rate  .  Then  system  reliability  Rs  can  be  written  as  a 
function  of  /„  t„  i  =  1 , 2,  •  •  • ,  k\  /.<?., 


^s(0  -  sty i>  4  4  4  >  J-h>  ri»  ?2*  •  •  •  >  lk) 

where  t,  =  /,(/)  is  the  operating  time  for  component  /,  /  =  1,2  ,k . 


Let  =  max{2„  A2,  ••• ,  >.*}  and  r,  =  /=  1, 2,  •  •  • ,  k.  Then  one  can  write 

UO  =  Si^-nv  rl )  *  4  *  >  rlo  hi  4  4  4  >  lk) 

A 

If  the  r,  were  known  and  were  an  upper  100(1  -  a)%  confidence  limit  for  ,  the 
corresponding  lower  confidence  limit  for  Rs(t)  would  be 

A  A 

Rs(!)l(z)  ~  SV-m,U(x)>  rh'f‘  j>'k>  hi  4  4  4  >  lk) 

Specifically,  if  the  system  is  a  series  system  of  independent  components,  so  that 


men, 


(2.1) 


If  n,  items  of  component  /  are  tested  until  failure,  T,  denotes  the  total  test  time  accumu- 

*  k 

lated  by  all  n,  items  and  n  =  then  2A„Y,r,T,  is^L*  See  Bain  and  Engelhardt  (Ref. 

IB]  1*5 1 

14]. 

An  upper  confidence  limit  for  is 
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(2.2) 


A 

^ m,U(a )  ~ 


where  xi»  is  the  100(1  -  a)th  percentile  of  the  x„  distribution.  Corresponding  equations 
for  truncated  testing  are  similar. 

If  the  r,  are  unknown,  the  following  methods  estimate  the  values  for  r,  from  the  data. 
One  method  uses  the.  likelihood  ratio  estimate  for  r,.  The  second  method  uses  a  jacknifed 

A 

version  of  the  first  method.  The  two  resulting  confidence  limits  Rsxw  and  Rs,lM  ,  with 
and  without  jacknifing  r,  respectively,  are  compared  for  relative  accuracy.  Quenouille 
[Ref.  15:  pp.  353-360]  first  reported  a  method  for  estimating  ratios  that  reduced  the  bias 
without  increasing  the  variance.  Miller  [Ref.  16:  pp.  1-15]  gives  an  excellent  review  of  the 
jacknife  method  which  includes  a  discussion  on  the  application  of  jacknifing  to  estimat¬ 
ing  ratios. 

A.  LOWER  CONFIDENCE  LIMIT  RSL  WITHOUT  JACKNIFING 

In  this  case,  the  maximum  likelihood  estimate  of  the  ratio  r,  =  -r1-  is 


A  A  A  A 

where  /.,-nJT,  and  =  max  {Xu  •  ••,/*).  The  resulting  approximate  upper  confr 
dcnce  limit  for  ).m  is 


where 


k 

n  =  ,  and 


(2.4) 


(2.5) 


and  T,j  denotes  the  failure  time  of  the  Jth  test  for  component  /. 
The  resulting  approximate  confidence  bound  Rs,u,y  is  given  by 
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(2.6) 


RS,L(*) 


B.  LOWER  CONFIDENCE  LIMIT  RStL  WITH  JACKNIFING 

A 

The  definitions  for  ,  r, ,  n, ,  and  T,  in  Section  A  are  also  used  in  this  section.  Let 

A 

/vW  denote  the  estimate  for  by  removing  Tt]  from  the  data  ;  /.<?., 


and 


2>„ 

/=! 


j  =  1,2,  •••,«;,  /#y 


(2.7) 


/W«  ~ 


j  =  1,2, • 


(2.8) 


Then  the  jacknifed  ratio  estimate  /•,*  is  given  by 


A* 


n  =  >/ 


{n  -  1)^ 
;=i 


where  «'  =  min(«,,  •  •♦,«*)  and 


A 

a 

}‘R')  0 


j  —  1,2,  •••,/! 


A 


mKr) 


Now  define  /J,.w  by 


v2 

A«,  2rt 

k 


£ 
i=  i 


The  corresponding  confidence  bound  is  given  by 


(2.9) 


(2.10) 


(2.H) 
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(2.12) 


C  k 

A  l  A*  \ 1  A* 

Rs,L(*)  =  exP)  “  aV{*)2~i  r‘ 

(.  i=\ 

When  the  nit  i  =  1, 2, •  k  differ  considerably,  this  jacknife  estimation  proce¬ 
dure  can  be  unbalanced.  That  is,  the  number  of  data  points  used  to  compute  the  jacknife 
estimate  will  differ  from  one  component  to  another.  It  was  decided  to  use  this  jacknifmg 
procedure  and  determine  the  effect  of  differing  sample  sizes  by  examining  the  results  of 
the  simulations. 

In  equation  (2.9),  the  jacknife  estimate  is  constructed  by  using  only  the  first  ri  ob¬ 
servations  from  each  component  to  obtain  rm  where  ri  =  minjrt,,  n2,  •  •  • ,  nk}.  Of 
course  it  is  rather  arbitrary  to  take  the  first  ri  observations  from  n, ,  since  any  ri  of  the 
collection  of  n,  values  could  be  used. 

C.  SIMULATIONS  AND  RESULTS  FOR  EXPONENTIAL  CASE 
1.  Simulation 

a.  Simulation  language  and  package 

The  programming  language  used  to  simulate  this  problem  was  VS 
FORTRAN  on  an  IBM  3033.  Also  LEXPN  in  LLRAXII  was  used  to  generate  observed 
exponential  random  variates.  SHSORT  was  used  to  perform  the  sort  routines. 

b.  Cases  and  Input  parameters 

The  six  input  parameters  below  determine  the  conditions  for  exercisin,.  *ne 
simulation  runs.  System  reliability  is  determined  by  five  parameters.  Values  of  these 
parameters  are  given  in  the  tables  that  show  the  results.  The  test  plan  simulated  was  to 
test  n,  items  until  all  fail. 

•  number  of  component  types,  k  :  5  and  15 

•  system  reliability,  Rs  :  0.9  and  0.975 

•  significance  level,  a  :  0.05  and  0.2 

•  component  time,  t,  :  small  to  large  (see  tables) 

•  sample  size,  :  small  to  large  (see  tables) 

k 

Reliability  of  a  series  system  is  expressed  as  Rs  —  exp{  -  £  /, We  chose  arbitrarily 

1=1 

to  determine  the  failure  rates,  X, ,  from  this  equation  by  assuming  all  X,t,  to  be  equal, 
consequently, 
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-lnflc 

Ai  =  ~nr  ■ 

c.  Replications 

The  procedure  was  replicated  1000  times  for  each  case  to  get  1000  values 

A 

of  Rs.lm  and  1000  values  of  RStLM  .  We  order  each  set  of  the  1000  values  of  RS<LM  and 

A 

R-s.lm  In  ascending  order.  Then  the  two  approximate  confidence  bounds,  i?s,iW{iooo(i-,))  and 

A 

R-s,lm  (iooc(i-,)> >  f°r  Rs  are  the  1000(1  -  oc)th  order  values  of  these  two  sets  of  data.  Finally, 
for  each  of  the  two  ordered  sets  of  data,  we  find  the  order  indices  jt  and  j2  for  which 

A 

Rs.n»uo  and  Rsm*xj2)  are  cl°sest  t0  Rs  for  their  respective  sets.  Then  jj  1000  and 
j2l  1000  are  called  the  two  corresponding  simulated  true  confidence  levels. 

2.  Results 

Tables  1  through  3  show  the  results  of  the  3  cases  simulated.  The  results  indi¬ 
cate  that  this  interval  estimation  method  using  estimates  of  failure  rate  ratios  will  yield 
quite  accurate  lower  confidence  limits  for  system  reliability  when  components  have  un¬ 
known  constant  failure  rates.  The  jacknife  method  also  yields  very  accurate  confidence 
limits  which  are  slightly  conservative.  That  is  the  100(1  —  a)  percentile  points  of  Rs,lM  , 
given  in  the  tables,  are  slightly  less  than  the  true  value  of  Rs.  Consquently, 
P(Rsm  <  Rs)  >  1  -  a  •  Alternatively,  one  can  say  RSMi)  is  a  conservative  100(1  -ot) 
percent  lower  confidence  limit  procedure  for  Rs  . 
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Table  1.  RELIABILITY  OF  A  SERIES  SYSTEM  WITH  SMALL  NUMBER 
_ (  LESS  THAN  10 )  OF  SAMPLE  SIZES  -  EXPONENTIAL  CASE 


Number  of 
Compo¬ 
nents 

Reliability 

Lower  Confidence 
Limit 

True  Confidence 
Limit 

of  System 

a 

■Eft'll 

WITH 

w/6 

WITH 

Jacknifing 

Jacknifing 

Jacknifing 

.90 

.2 

0:9017 

0.8904 

0,7720 

0.9100 

5 

.05 

0.8985 

0.8889 

0.9620 

0.9780 

.975 

.2 

0.9750 

0.9720 

0.8020 

.05 

0.9746 

0.9729 

0.9600 

0.9800 

.90 

.2 

0.9019 

0.8901 

0.7600 

0.9500 

15 

.05 

0.9012 

0.8933 

0.9300 

0.9800 

.975 

.2 

0.9745 

0.9724 

0.8300 

0.9600 

.05 

0.9746 

0.9725 

0.9600 

0.9800 

Sample  sizes  for  5  components  are  9,  7,  10,  8,  6 

Sample  sizes  for  15  components  are  6,  7,  5,  6,  7,  8,  9,  5,  8,  7,  9,  10,  7,  9,  6 
Component  times  t,  are  2,  3,  7,  8,  10,  3,  7,  10,  1,  7,  8,  3,  10,  1,  8 


Table  2.  RELIABILITY  OF  A  SERIES  SYSTEM  WITH  MEDIUM  NUMBER 
(  LESS  THAN  30  )  OF  SAMPLE  SIZES  -  EXPONENTIAL  CASE 


Number  of 
Compo¬ 
nents 

Reliability 

Lowe;  Confidence 
Limit 

True  Confidence 
Limit 

of  System 

oc 

W;0 

Jacknifing 

WITH 

Jacknifing 

W;0 

Jacknifing 

WITH 

Jacknifing 

.90 

.2 

0.8999 

0.8978 

0.8020 

0.8660 

5 

.05 

0.8990 

0.8983 

0.9600 

0.9680 

.975 

.2 

0.9750 

0.9745 

0.8000 

0.8460 

.05 

0.9748 

0.9745 

0.9620 

0.9740 

.90 

.2 

0.8992 

0.8960 

0.8300 

0.9100 

15 

.05 

0.8980 

0.8981 

0.9600 

0.9700 

.975 

.2 

0.9751 

0.9738 

0.7900 

0.9000 

.05 

0.9750 

0.9741 

0.9500 

0.9800 

Sample  sizes  for  5  components  are  30,  21,  10,  15,  26 

Sample  sizes  for  15  components  are  3,  7,  10,  15,  20,  15,  7,  5,  30,  20,  10,  7,  13,  21,  30 


Table  3.  RELIABILITY  OF  A  SERIES  SYSTEM  WITH  LARGE  NUMBER 


(  LESS  THAN  100  )  OF  SAMFLE  SIZES  -  EXPONENTIAL  CASE 


Number  of 
Compo¬ 
nents 

Reliability 
of  System 

a 

Lower  Confidence 
Limit 

True  Confidence 
Limit 

W/O 

Jacknifing 

WITH 

Jacknifing 

W/O 

Jacknifing 

WITH 

Jacknifing 

5 

.90 

.2 

0.8999 

0.8991 

0.8080 

0.8200 

.05 

0.8990 

0.8992 

0.9680 

0.9620 

.975 

.2 

0.9750 

0.9750 

0.8540 

0.8060 

.05 

0.9751 

0.9750 

0.9380 

0.9520 

15 

.90 

.2 

0.9000 

0.8988 

0.8010 

0.8750 

.05 

0.9000 

0.S996 

0.9500 

0.9590 

.975 

.2 

0.9750 

0.9747 

0.8460 

0.8730 

.05 

0.9750 

0.9750 

0.9520 

0.9470 

Sample  sizes  for  5  components  are  30.  63,  75,  98,  26 

Sample  sizes  for  15  components  are  15,  40,  35,  17,  26,  67,  50.  65,  80,  32,  95,  100,  15, 
45,  30 


III.  INTERVAL  ESTIMATION  PROCEDURE  -  WEIBULL  CASE 

A.  LOWER  CONFIDENCE  LIMIT  RyL 

Consider  a  series  system  with  k  components.  Let  the  time  to  failure  X,  of  compo¬ 
nent  /  have  a  Weibull  distribution  with  density 

m  -  >!< v <!<-' t* p{ ,  r,>o.  p.i) 

Then 

=  exp{-(V/']  ,  (f>0  (3.2) 

and 


Rs{!)  =  exp 


t>0  .  (3.3) 


where  /'  —  /A ,  ).'n  -  max/.;  and  r,  =  If  the  /?,  are  known,  X?>  has  constant 

failure  rate  /.?'  and  the  procedures  in  Chapter  II  can  be  used  to  obtain  Rs,lw  anc* 
Rs.  LM  with  T,j  replaced  by  7^  in  equation  (2.5).  It  is  well  known  that  Y  =  In  X,  has 
an  extreme  value  distribution  with  CDF 


F}{v)  -  1  -  exp  j  — 


where  =  ln(l//.,)  and  <5,  =  1  //?,. 

Engelhardt  and  Bain  [Ref.  17:  p.  323]  have  developed  the  following  simple  unbiased 
estimators  for  c,  and  5, ,  using  ordered  values  YtJ  =  In  XI(J) 


(3.4) 


where  s  —  [0.84/7,]  s  largest  integer  <  0.84  n,  and  is  the  jth  order  statistics  from  the 
sample  of  size  n,  of  X, .  Also 
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(3.6) 


where  y  =  0.5772  and  yt-YJY,I!nt .  Let 


and 

A 

Ttj  -  Xfj‘  /  =»  1,2, •••,«/  y  -  1,2,... ,A.  (3.7) 

We  approximate  the  distribution  of  TtJ  by  the  exponential  distribution  with  failure  rate 

A 

Xft  =  ?.]  and  proceed  as  in  Chapter  II.  Define 


'•/ 


!1l 

T, 


"I  A  A 

where  T,  =  £ T,.  i  =  1,2, Let  =  max /,*  and 

.=i  < 


Then  an  approximate  upper  confidence  limit  for  is  given  by 


o* 


(3.8) 


(3.9) 


(3.10) 


and  the  corresponding  approximate  lower  confidence  limit  RSM  for  Rs(t)  is  given  by 


(3.11) 


k 

where  n  =  . 

i=i 

This  procedure  is  labeled  the  Formula  procedure  in  the  tables  that  follow  in  this 

A 

section.  Its  distinguishing  feature  is  the  equation  for  /?,  given  by  equation  (3.6).  An  al¬ 
ternative  procedure,  labeled  the  Newton  -  Raphson  procedure,  estimates  /?,  using  the 


maximum  likelihood  procedure  which  is  obtained  using  a  Newton  -  Raphson  approxi- 

A 

mation  method.  Equations  for  ft  under  the  Newton  -  Raphson  procedure  are  provided 
in  Appendix  A. 

B.  SIMULATIONS  AND  RESULTS 
1.  Simulation 

a .  Simulation  language  and  package 

The  programming  language  used  to  simulate  this  problem  was  VS 
FORTRAN  on  an  IBM  3033.  Also  LEXPN  in  LLRANII  was  used  to  generate  observed 
exponential  random  variates.  SHSORT  wras  used  to  perform  the  sort  routines. 

b.  Cases  and  Input  parameters 

The  six  input  parameters  below  determine  the  conditions  for  exercising  the 
simulation  runs.  System  reliability  is  determined  by  five  parameters.  Values  of  these 
parameters  are  given  in  the  tables  that  show  the  results.  The  test  plan  simulated  w’as  to 
test  n,  items  until  all  fail. 

•  number  of  component  types,  k  :  5  and  15 

•  system  reliability,  Rs  :  0.9  and  0.975 

•  significance  level,  a  :  0.05  and  0.2 

•  component  time,  r,  :  small  to  large  (see  tables) 

•  sample  size,  n,  :  small  to  large  (see  tables) 

k  _ 

Reliability  of  a  series  system  is  expressed  as  Rs  =  exp{  -  £  •  Failure  rates  can 

i=  1 

be  determined  from  that  equation  by  assuming  all  X,i,  to  be  equal,  consequently, 

,  _  ~ln  Rs 
k  t, 

c.  Replications  ■ 

The  procedure  was  replicated  1000  times  for  each  case  to  get  1000  values 

A 

of  Rs<Ui)  and  1000  values  of  RSiLM  .  We  order  each  set  of  the  1000  values  of  RSiLM  and 

A 

Rs,u>)  ascending  order.  Then  the  twro  approximate  confidence  bounds,  R$,i(,)<iooo(i-.»  and 
Rs  ars  the  1000(1  -  a)th  order  values  of  these  two  sets  of  data.  Finally, 
for  each  of  the  two  ordered  sets  of  data,  we  find  the  order  indices  j\  and  j2  for  wfhich 

A 

Rs.imv j)  a°d  Rs.lw 2)  are  the  cl°sest  to  Rs  for  their  respective  sets.  Then  jj  1000  and 
j2 1  1000  are  called  the  two  corresponding  true  confidence  levels. 
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2.  Results 

Tables  4  and  5  display  the  results  of  the  simulations  and  determine  the  accuracy  of 

A 

RSiLM  §ivcn  m  equation  (3. 1 1)  as  a  lower  confidence  limit  procedure  for  system  reliability, 
Rs,  for  parameter  values  /?„  t ,  and  Rs  given  in  the  tables.  The  terms  in  the  tables  have 
the  same  meaning  as  the  corresponding  terms  in  Tables  1  through  3  which  were  ex¬ 
plained  in  Section  2.C.2  .  The  procedure  would  be  exact  for  the  Formula  method  if  the 
values  in  the  Formula  column  equal  the  corresponding  numbers  in  the  same  row  in  the 
Reliability  of  System  column.  In  Table  4  for  example,  the  80th  percentile  point  of  .9205 
for  the  Newton  -  Raphson  procedure  is  more  accurate  than  the  Formula  procedure 
which  has  an  80th  percentile  point  of  .9324.  The  last  column  shows,  however,  that  what 
we  have  called  an  80%  lower  confidence  limit  procedure  is  in  fact  closer  to  a  61%  pro¬ 
cedure.  The  accuracy  improves  if  the  sample  size  is  increased  from  15  to  30  as  indicated 
in  Table  5.  The  accuracy  improves  even  more  if  the  number  of  compom.  :ts  in  the  system 
increases  from  5  to  15  as  indicated  in  Tables  4  and  5. 

A 

The  results  in  Tables  4  and  5  show  that  the  RUl)  method  given  by  equation 
(3.11)  is  too  optimistic  -  especially  for  small  sample  sizes  and  for  systems  with  a  small 
number  of  components.  Modifications  to  this  procedure  are  needed.  These  tables  also 
indicate  that  the  Newton  -  Raphson  method  is  superior  to  the  Formula  method. 
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Table  4.  ACCURACY  OF  Rsl  AS  A  LOWER  CONFIDENCE  LIMIT  FOR  Rs 
WITH  SAMPLE  SIZES  OF  15  -  WEIBULL  CASE 


Number  of 
Compo¬ 
nents 

Reliability 
of  System 

a 

Lower  Confidence 
Limit 

True  Confidence 
Limit 

Formula 

New  on 
Raphson 

Formula 

Newton 

Raphson 

5 

.90 

.2 

0.9324 

0.9205 

0.4820 

0.6060 

.05 

0.9512 

0.9406 

0.5380 

.975 

.2 

0.9849 

0.9824 

0.5040 

0.6540 

.05 

0.9905 

0.9S91 

0.5600 

0.6660 

15 

.90 

.2 

0.9193 

0.9030 

0.5440 

0.7620 

.05 

0.9416 

0.9288 

0.6060 

0.7760 

.975 

.2 

0.9802 

0.9747 

0.6160 

0.8080 

.05 

0.9865 

0.9S28 

0.6580 

0.8240 

Initial  beta  for  5  components  are  2.2,  2.4,  2.6,  2.8,  2.5 

Initial  beta  for  15  components  are  2.2,  2.4,  2.6,  2.S,  2.5,  2.3,  2.7,  2.5,  2.2,  2.6,  2.9,  2.8, 
2.4,  2.6,  2.1 

Component  times  t,  are  2,  3,  7,  8,  10,  3,  7,  10,  1,  7,  8,  3,  10,  1,  8 


Table  5.  ACCURACY  OF  R5L  AS  A  LOWER  CONFIDENCE  LIMIT  FOR  Rs 
WITH  SAMPLE  SIZES  OF  30  -  WEIBULL  CASE 


Number  of 
Compo¬ 
nents 

Reliability 
of  System 

a 

Lower  Confidence 
Limit 

True  Confidence 
Limit 

Formula 

Newton 

Raphson 

Formula 

Newton 

Raphson 

5 

.90 

.2 

0.9249 

0.9125 

0.4100 

0.6320 

.05 

0.9424 

0.9321 

0.5100 

0.6940 

.975 

.2 

0.9837 

0.9798 

0.4380 

0.6160 

.05 

0.9S79 

0.9845 

0.4600 

0.6760 

15 

.90 

.2 

0.9207 

0.9036 

0.5560 

0.7200 

.05 

0.9333 

0.9198 

0.4380 

0.7520 

.975 

.2 

0.9808 

0.9758 

0.5120 

0.7520 

.05 

0.9850 

0.9810 

0.5200 

0.7920 

Initial  beta  for  5  components  are  2.2,,  2.4,  2.6,  2.8,  2.5 

Initial  beta  for  15  components  are  2.2,  2.4,  2.6,  2.8,  2.5,  2.3,  2.7,  2.5,  2.2,  2.6,  2.9,  2.8, 
2.4,  2.6,  2.1 


14 


IV.  INTERVAL  ESTIMATION  PROCEDURE  -  NORMAL  CASE 


A.  BACKGROUND 

In  recent  years  it  has  become  popular  to  model  mechanical  reliability  as  P(X  >  Y)  where 
X  denotes  strength  and  Y  denotes  stress.  This  formulation  of  reliability  is  an  average  of 
P(A'>>-|  Y=y) ,  since 

P(X>Y)=EyP(X>y)  =  jp(X>y)fy(y)dy  , 

An  alternative  model  is  the  worst  case  approach.  In  the  worst  case  model,  one  se¬ 
lects  a  worst  case  value  ofy,  sayy0  ,  and  then  designs  the  strength,  X , of  the  component 
so  that  P(X  >y0 )  =  R0  where  R0  is  a  reliability  requirement.  If  X  has  a  normal  distrib¬ 
ution  then  this  requirement  imposes  constraints  on  the  mean  and  variance  of  X  .  This 
is  usually  done  in  a  manner  to  comply  with  standard  "safety  factor"  procedures. 

The  average  model,  P(X  >  Y),  uses  two  random  variables  and  is  subject  to  more 
random  error  than  the  worst  case  model.  The  accuracy  of  the  expression  P{X  >  Y)  for 
values  of  this  expression  close  to  one  is  highly  suspect  due  to  deviations  in  the  tail 
probabilities  of  both  X  and  Y  from  those  assumed  in  the  model.  It  is  common  to  assume 
that  both  X  and  Y  have  normal  distributions.  Truncated  normal  distributions  would  be 
more  appropriate  for  many  types  of  mechanical  equipment.  Harris  and  Soms  [Ref.  18: 
pp.  650-663}  discuss  implications  of  this  problem.  Very  significant  errors  in  point  and 
interval  estimation  for  reliability  are  readily  demonstrated  when  X  is  truncated  normal 
but  assumed  to  be  normal  in  the  more  simple  model  which  specifies  that  R  =  P(X  >  Y). 
Table  10  shows  this  effect  when  X  is  truncated  above  at  n  +  1.645<7 ,  where  Z,  is  the  100 
(1  -  ec)th  percentile  point  of  the  standard  normal  distribution.  Church  and  Harris  [Ref. 
12:  pp.  49-54]  and  Downton  [Ref.  13:  pp.  551-558]  have  developed  approximate  confi¬ 
dence  intervals  for  P( X  >  Y)  when  X  is  normal  with  unknown  mean  and  variance  and  Y 
has  the  standard  normal  distribution. 
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B.  EXACT  AND  APPROXIMATE  INTERVAL  ESTIMATES 

Minimum  variance  unbiased  estimators  (MVUE)  for  R  =  P(X>y)  are  well  known 
when  X  is  normal  with  unknown  mean  and  known  variance  and  also  when  the  variance 
is  unknown.  In  the  former  case,  Lieberman  and  Resnikoff[Ref.  8]  developed  the  result 


A  ,  X  -  v 
R  =  d> 


which  is  MVUE  for  R  where  is  the  standard  normal  cumulative  distribution  function. 
When  the  variance  is  unknown,  several  versions  of  an  integral  expression  for  P(X>y ) 
have  been  developed  by  Lieberman  and  ResnikofT[Ref.  8] ,  Basu  [Ref.  19]  and  Folks  and 
others  [Ref.  9]. 

Exact  lower  confidence  interval  estimates  for  P(X>y)  when  X  is  normal  with  un¬ 
known  mean  and  variance  involve  the  non-central  t  -  distribution.  Owen  and  Hua  [Ref. 
11]  developed  tables  of  the  lower  90%  and  95%  confidence  limit  values  RL  for  P{X>y) 
based  on  the  non-central  /  -  distribution.  These  values  are  tabled  for  values  of  k  in  the 
range  -3.0  (.2)  6.0  and  sample  sizes  n  =  2  (1)  18,  21  (3)  30,  40  (20)  100  ,  where 
k  =  (.v  -  y)js  and  x  and  s  are  the  sample  mean  and  standard  deviation.  Their  tables  are 
essentially  exact.  An  approximation  to  their  exact  tabled  values  is  given  by  R'L  where 
Rl  =  <!>(’/) , 


V  =  k  — 


■  +  ■ 


2{n  -  Jk ) 


Z3,  fl-1 


(4.1) 


k  =  (.v  —y)/s  and  r>in_,  is  the  100(1  -  a)th  percentile  of  the  t  distribution  with  n  —  1  de¬ 
grees  of  freedom.  Equation  (4.1)  was  developed  in  this  thesis.  It  is  an  extensive  ad  hoc 
modification  of  an  equation  developed  by  Church  and  Harris  [Ref.  12].  Tables  6  and  7 
display  values  of  RL  =  <h(y) ,  y ,  /,  /?].  ,  and  R[  -  RL  for  k  -  1,  2,  3,  4,  sample  sizes  n 
=  10,  IS,  30  and  confidence  levels  90%  and  95%  .  RL  and  y  denote  the  "exact"  lower 
confidence  limits  and  corresponding  0_l(RJ  values  from  Owen  and  Hua  [Ref.  11].  Both 
y*  and  R'L  are  given  by  equation  (4.1).  The  accuracy  of  the  approximate  confidence  in¬ 
terval  is  quite  good  relative  to  the  values  for  RL  given  by  Owen  and  Hua  [Ref.  1 1]. 
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Table  6.  APPROXIMATE  (  Rl  )  AND  EXACT  (  RL  )  90%  CONFIDENCE 
_ LIMITS  FOR  P(X >  Y0)  _ 


n 

m 

mm 

RiMy) 

V 

y* 

*[:<!>('/) 

M ebbb 

10 

1.3830 

l 

mamm 

0.47194 

0.45454 

0.67528 

-0.00628 

X. 

0.89130 

1,23397 

1.20199 

0.88517 

-0.00613 

3 

0.97453 

1.95262 

1.88991 

0.97025 

-0.00428 

4 

0.99602 

2.65302 

2.54950 

0.99437 

-0.00165 

18 

1.3334 

1 

0.73037 

0.61385 

0.61133 

0.72950 

-0.00087 

2 

0.92569 

1.44512 

1.44038 

0.92488 

-0.00081 

3 

0.98755 

2.24314 

2.23150 

0.9Su84 

-0.00071 

4 

0.99877 

3.02679 

3.00614 

0.99859 

-0.00018 

30 

1.3114 

1 

0.75937 

0.70424 

0.70508 

0.75960 

0.00023 

2 

0.94256 

1.57740 

1.57852 

U. 94248 

-0.00008 

3 

0.99233 

2.42404 

2.42459 

0.99205 

-0.00028 

mm 

0.99945 

3.26207 

3.25926 

0.99940 

-0.00005 

Table  7.  APPROXIMATE  (  Rl  )  AND  EXACT  (  RL  )  95%  CONFIDENCE 
_ LIMITS  FOR  P(X>  Y0) _ _ 


n 

Ao.'.r-l 

99 

y 

y‘ 

/&<&(/) 

ma 

10 

1.8331 

B9 

0.63052 

0.33311 

0.27702 

0.60912 

-0.02140 

2 

0.85187 

1.04477 

0.94228 

0.82692 

-0.02495 

3 

0.95565 

1.70308 

1.52864 

0.93655 

-0.01910 

mm 

0.99031 

2.33813 

2.07743 

0.98075 

-0.00956 

18 

1.7396 

i 

0.69504 

0.51007 

0.49292 

0.68896 

-0.00608 

2 

0.9034S 

1.30221 

1.26991 

0.89777 

-0.00571 

3 

0.97996 

2.05344 

1.99739 

0.97674 

-0.00322 

KB 

0.99735 

2.78717 

2.70338 

0.99638 

-0.00097 

30 

1.6991 

1 

0.7336! 

0.62369 

0,61789 

0.73167 

-0.00194 

2 

0.92855 

1.46579 

1.45391 

0.92677 

-0.00178 

3 

098855 

2,27522 

2.25448 

0.98758 

-0.00097 

4 

0.99894 

3.07140 

3.04028 

0.99873 

-0.00021 
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Tables  8  and  9  display  the  results  of  computer  simulations  with  1,000  replications 
to  check  the  accuracy  of  the  R'L  method  for  80%  and  90%  lower  confidence  limits  for 
P(X>y)  =  R  fory  =  3,  with  various  values  of  a  and  n  determined  so  that  R  equals 
the  values  shown.  The  procedure  would  be  exact  at  the  80%  level  if  the  values  in  the 
column  labelled  a  =  .2  equal  the  corresponding  values  of  R  in  the  same  row.  The  "true" 
confidence  level  corresponds  to  the  index  i(R)  of  1,000  ordered  values  of  R'L{,}  for  which 
R'LWA R)  =  R  .  For  example,  in  the  seventh  row  of  Table  8,  R  =  .950  ,  ax  -  20 ,  N  =  10, 
Rlmm  was  -9575,  and  R'Lmm  was  .9533.  Also  Rl(.20)rSt  =  -950  and  R'Lmm  =  .950. 
Tables  8  and  9  indicate  that  the  R'L  procedure  given  by  equation  (4.1)  is  quite  accurate 
at  the  80%  level  of  confidence  for  the  cases  simulated. 
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Table  8.  ACCURACY  ANALYSIS  OF  R*L  PROCEDURE  FOR  80%  AND  90% 
LOWER  CONFIDENCE  LIMITS  FOR  P{X>  3)  WHEN  X  IS 
NORMALLY  DISTRIBUTED 


R 

n 

^Z..1000(l->> 

True  Confidence  Level  J 

a  =  .2 

a  =  .1 

a  =  .2 

a  =  .1 

10 

0.9624 

0.9504 

0.7310 

0.8970 

0.5 

25 

0.9546 

0.9542 

0.7560 

0.8820 

75 

0.9517 

0.9501 

0.7S00 

0.8980 

10 

0.9606 

0.9531 

0.7470 

.950 

1.0 

25 

0.9501 

0.9552 

0.7970 

0.8770 

75 

0.9509 

0.9528 

0.7900 

0.8770 

10 

0.9575 

0.9533 

0.7580 

0.8860 

20.0 

25 

0.9558 

0.9551 

0.7410 

0.8730 

75 

0.9511 

0.9511 

0.7830 

0.8820 

10 

0.9924 

0.9913 

1.0000 

0.8840 

0.5 

25 

0.9906 

0.9898 

0.7S00 

0.9070 

75 

0.9903 

0.9907 

0.9800 

0.8780 

10 

0.9921 

0.9S99 

0.9990 

0.9000 

.990 

1.0 

25 

0.9904 

0.7570 

0.8910 

75 

0.9904 

0.9898 

0.7S00 

0.9060 

10 

0.9923 

0.9922 

0.9880 

0.8750 

20.0 

25 

0.9909 

0.9908 

0.7760 

0.8800 

0.9903 

0.7940 

0.8920 

10 

0.9966 

0.994S 

1.0000 

0.9030 

0.5 

25 

0.9955 

0.9949 

0.7790 

0.9020 

75 

0.9956 

0.9949 

0.7540 

0.9060 

10 

0.9966 

0.9952 

1.0000 

0.S950 

.995 

1.0 

25 

0.9960 

0.9955 

0.8630 

0.8880 

75 

0.9949 

0.9947 

0.8040 

0.9150 

10 

0.9970 

0.9951 

1.0000 

1.0000 

20.0 

25 

0.9960 

0.9957 

0.7720 

0.8700 

75 

0.9951 

0.9948 

0.7940 

0.9140 
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Table  9.  ACCURACY  ANALYSIS  OF  R*L  PROCEDURE  FOR  80%  AND  90% 
LOWER  CONFIDENCE  LIMITS  FOR  P(X>30)  WHEN  X  IS 
NORMALLY  DISTRIBUTED 


R 

WM 

^£.1000(1 -s> 

True  Confidence  Level 

mm 

n 

<x  =  .2 

a  =  .1 

a  =  .2 

a  =  .1 

10 

0.9621 

0.9466 

0.7510 

0.9000 

0.5 

25 

0.9530 

0.9519 

0.7680 

0.8890 

75 

0.9507 

0.9491 

0.7880 

0.9070 

10 

0.9605 

0.9528 

0.7480 

0.8910 

.950 

1.0 

25 

0.9497 

0.9549 

0.8030 

0.8810 

75 

0.9509 

0.9527 

0.7910 

0.8780 

10 

0.9575 

0.9533 

0.7580 

0.8860 

20.0 

25 

0.9558 

0.9551 

0.7410 

0.8730 

75 

0.9511 

0.9511 

0.7830 

0.8820 

10 

0.9923 

0.9911 

1.0000 

0.8860 

0.5 

25 

0.9900 

0.9890 

0.7990 

0.9140 

75 

0.9901 

0.9904 

0.79S0 

0.8860 

10 

0.9921 

0.9899 

0.9990 

0.9310 

.990 

1.0 

25 

0.9911 

0.9902 

0.7620 

0.8950 

75 

0.9904 

0.9898 

0.7S00 

0.9070 

10 

0.9923 

0.9922 

0.98S0 

0.8750 

20.0 

25 

0.9909 

0.9908 

0.7760 

O.S800 

75 

0.9902 

0.9903 

0.7940 

0.8840 

10 

0.9965 

0.9946 

1.0000 

0.9040 

0.5 

25 

0.9952 

0.9946 

0.7900 

0.9140 

75 

0.9955 

0.9948 

0.7620 

0.9100 

10 

0.9966 

0.9951 

1.0000 

0.9280 

.995 

1.0 

25 

0.9959 

0.9953 

0.7610 

0.8920 

75 

0.9949 

0.9947 

0.8040 

0.9140 

10 

0.9970 

0.9951 

1.0000 

— 1 

20.0 

25 

0.9960 

0.9957 

0.7600 

0.8700 

75 

0.9951 

0.9948 

0.7940 

0.9340 
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Table  10  displays  the  inaccuracies  of  the  R’L  lower  confidence  limit  procedure  for 
P(X>  3)  using  equation  (4.1),  which  assumes  A'  is  N(ji,  o2),  when  in  fact  X  has  the  dis¬ 
tribution  of  a  normal  random  variable  with  mean  n  and  variance  o 2  that  has  been  trun¬ 
cated  at  n  +  1.645cr.  Note  that  /;  +  1.6450  >  >  3,  because  P(X  >  3)  =  .95,  .99  and  .995. 

Examination  of  Table  10  reveals  gross  inaccuracies;  consequently,  even  when  the 
distribution  of  X  is  truncated  far  above  the  value  y,  the  exact  lower  confidence  limit  for 
P( X  >y)  can  be  greatly  in  error  when  computed  under  the  assumption  that  X  is  normal. 
This  problem  will  be  compounded  when  one  is  computing  "exact"  confidence  intervals 
for  P{X  >  Y)  assuming  both  X  and  Y  are  normal  when  in  fact  one  or  both  may  be  trun¬ 
cated  normal.  This  suggests  that  modeling  mechanical  reliability  as  P(X  >  y)  may  be 
more  risky  than  more  standard  engineering  app^oachs  for  modeling  mechanical  reliabil¬ 
ity  which  include  the  notion  of  safety  factors. 
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Table  10.  ACCURACY  ANALYSIS  OF  R*L  PROCEDURE  FOR  80%  AND  90% 
LOWER  CONFIDENCE  LIMITS  FOR  P(.Y>  3)  WHEN  X  IS  TRUN¬ 
CATED  NORMAL 
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V.  CONCLUSIONS  AND  RECOMMENDATIONS 


The  lower  interval  estimation  procedure  for  reliability  of  coherent  systems  which 
was  developed  in  Chapter  II  appears  to  be  accurate,  easy  to  use,  and  applicable  to  co¬ 
herent  systems.  Although  this  procedure  assumes  that  failure  times  of  all  components 
of  the  system  are  independent  and  have  exponential  probability  distributions,  it  can  be 
easily  extended  to  systems  with  component  failure  distributions  that  can  be  transformed 
into  exponential  failure  distributions  ;  e.g.,  Weibull  distribution  with  known  shape  pa¬ 
rameter.  This  procedure  has  potential  for  being  combined  with  a  similar  procedure  for 
systems  with  cyclical  components.  The  combined  procedure  would  provide  for  use  of 
binomial  component  test  data  and  exponential  component  test  data  to  compute  lower 
confidence  limits  on  the  reliability  of  coherent  systems  with  both  cyclic  and  continuous 
components.  Such  an  extension  could  use  a  failure  rate  ratio  estimation  procedure  sim¬ 
ilar  to  that  developed  in  this  thesis.  Such  a  method  should  be  explored. 

The  interval  estimation  method  for  the  reliability  of  a  system  with  components  that 
have  Weibull  failure  times  is  not  sufficiently  accurate  to  be  applied  to  systems  that  have 
10  or  fewer  components  each  with  ten  or  fewer  tests.  This  procedure  needs  further  study 
and  refinement. 

The  approximate  lower  confidence  limit  for  component  reliability  P(X>y)  when  X 
is  normally  distributed  with  unknown  mean  and  variance  is  very  accurate.  It  does  not 
require  an  extensive  set  of  tables  such  as  those  developed  by  Owen  and  Hua[Ref.  1 1], 
but  only  requires  the  use  of  the  standard  /  tables. 
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APPENDIX  A.  MLE  OF  WE1BULL  PARAMETERS  BY  THE  NEWTON  - 

RAPHSON  METHOD 

Let  X  ~  lVEI{a,  /?).  Then  the  likelihood  function  for  the  first  r  ordered  observations 
from  a  random  sample  of  size  n  is  given  by 


r 


Setting  the  partial  derivatives  of  this  likelihood  with  respect  to  a  and  /?  equal  to  zero 
gives  the  MLE's  a  and  /?  as  solutions  to  the  equations 


and 


A 


a  = 


where  n  =  r.  It  can  be  shown  that  these  equations  have  unique  solutions  which  are  the 
maximum  likelihood  estimates.  The  NEWTON  -  RAPHSON  method  for  solving  an 

A  A 

equation  g(P)  =  0  is  to  determine  successive  approximations  /?,,  where 

A  A  A  A 

=  fij  —  g(pj)lg'(pj)-  Therefore,  the  estimates  of  a  and  /?  can  be  solved  by  letting 
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A 

The  derivative  of  g(/?)  is 


m  = 


4  ln*,)J 
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APPENDIX  B.  FORTRAN  CODE  FOR  INTERVAL  ESTIMATION 
PROCEDURE  ■  EXPONENTIAL  CASE 

PROGRAM  EXPONE 
C 

Q  ****************************************************************** 

c  *  * 

C  *  THIS  IS  A  PROGRAM  TO  COMPUTE  THE  TRUE  CONFIDENCE  LIMIT  OF  * 

C  *  A  SERIES  SYSTEM  WITH  /  WITHOUT  JACKNIFE  METHOD  AND  * 

C  *  COMPARE  THE  DIFFERENCE  OF  THOSE  RESULTS.  * 

C  *  * 

C  *  BELOWS  ARE  GIVEN  OR  ASSUMED.  * 


c 

* 

* 

c 

* 

N 

NUMBER  OF  COMPONENTS  TYPE 

* 

c 

* 

RSYS 

RELIABILITY  OF  A  SERIES  SYSTEM 

* 

c 

* 

ALPHA 

SIGNIFICANCE  LEVEL 

* 

c 

* 

TIME 

TEST  TIME 

* 

c 

c 

* 

* 

SAMPLE 

SAMPLE  SIZE  FOR  EACH  TYPE  OF  COMPONENT 

* 

* 

c 

c 

r i 

it 

* 

* 

* 

THESE  ARE  VARIABLES  USED. 

* 

* 

* 

0 

c 

B 

TEMPORARY  ARRAY  FOR  EXPONENTIAL  RANDOM  VARIATE 

* 

c 

B I GLAM 

LARGEST  VALUE  OF  LAMBDA 

* 

c 

it 

CASE 

COUNTER 

* 

c 

it 

CHISQR 

CHI-SQUARE  VALUE  FOR  GIVEN  ALPHA  AND  DF 

* 

c 

it 

CLR1 

TRUE  CONFIDENCE  LIMIT  FROM  NON-JACKNIFING 

* 

c 

it 

CLR2 

TRUE  CONFIDENCE  LIMIT  FROM  JACKNIFING 

* 

c 

* 

CTIME 

TEST  TIME  OF  INDIVIDUAL  COMPONENT 

* 

c 

* 

DF 

DEGREE  OF  FREEDOM  FOR  CHI-SQUARE. 

* 

c 

* 

DIFF 

ABSOLUTE  VALUE  OF  DIFFERENCE  BETWEEN  R1  &  R2 

* 

c 

* 

LAMBDA 

FAILURE  RATE  OF  EACH  COMPONENT  TYPE 

* 

c 

it 

LAMHAT 

LAMBDA  HAT  FOR  JACKNIFING 

* 

c 

it 

LAMHM 

FINAL  LAMBDA  HAT  FROM  NON-JACKNIFING 

* 

c 

it 

LAMHST 

FINAL  LAMBDA  HAT  FROM  JACKNIFING 

* 

c 

it 

LAMMAX 

LARGEST  VALUE  OF  LAMHAT 

* 

c 

it 

LOMIT 

LAMBDA  WITH  1  COMPONENT  OMITTED 

* 

r 

it 

KEY1 

TEMPORARY  ARRAY 

* 

c 

it 

KEY2 

TEMPORARY  ARRAY 

* 

c 

it 

R 

INITIALLY  COMPUTED  RELIABILITY  FOR  EACH 

* 

c 

it 

COMPONENT  TYPE 

* 

c 

it 

R1 

COMPUTED  RELIABILITY  BY  NON-JACKNIFING 

* 

c 

it 

R2 

COMPUTED  RELIABILITY  BY  JACKNIFING 

* 

c 

it 

RATIO 1 

RATIO  OF  LAMBDA  FOR  NON-JACKNIFING 

* 

c 

it 

RATI02 

RATIO  OF  LAMBDA  FOR  JACKNIFING 

* 

c 

it 

ROMIT 

RATIO  WITH  1  COMPONENT  OMITTED 

* 

c 

it 

ROMSUM 

USED  FOR  JACKNIFE,  SUM  OF  R  WITH  OMIT  1 

* 

c 

it 

RSTAR 

FINALLY  COMPUTED  R  VALUE  (  R  STAR  ) 

* 

c 

it 

RVAL1 

R(  500  *  (1-ALPHA)  )  FOR  NON-JACKNIFING 

* 

c 

it 

RVAL2 

R(  500  *  (1-ALPHA)  )  FOR  JACKNIFING 

* 

c 

it 

T 

TOTAL  TEST  TIME  OF  EACH  COMPONENT  TYPE 

* 

c 

it 

ZALPHA 

RIGHT  PERCENTILE  POINTCNORMAL  DISTRIBUTION) 

* 
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C 

C 

PARAMETER  (NN  =  500) 

REAL  TIME(15),  ALPHA(2) ,  ZALPHA(2),  RSYS(2),  LAMBDA( 15) ,  P 
REAL  LAMMAX,  LAMHAT( 15),  B(100),  RSTAR(15),  LAMHM 
REAL  RATI01C15),  RATI02(15),  BIGLAM,  DIFF(15) 

REAL  LOMIT( 15 , 100) ,  CTIME( 15 , 100) ,  ROMIT( 15 , 100) ,  T( 15) 

REAL  ROMSUM  ,  SUM1,  SUM2,  LAMHST,  Rl(NN) ,  R2(NN) ,  KEYl(NN) 

REAL  RVALK15),  RVAL2(15),  CLR1(15),  CLR2(15),  KEY2(NN) 
REAL  DIFFRK 15) ,  DIFFR2(15),  RC 15) 

INTEGER  SAMPLE (15),  N(2),  DF,  CASE 

DATA  N  /  5,0  / 

DATA  RSYS  /  .9,  .975  / 

DATA  TIME  /  2,3,7,8,10,3,7,10,1,7,8,3,10,1,8  / 

DATA  SAMPLE  /  30,63,75,98,26,15,7,5,30,20,10,7,13,21,30  / 
DATA  ALPHA  /  .2,  .05  / 

DATA  ISEED  /  1736  / 


C  /*  SUBROUTINE  ZTABLE  COMPUTES  THE  RIGHT  PERCENT  POINT  ZALPHA  */ 

C  /*  FROM  RIGHT  CUMULATIVE  PROBABILITY  ALPHA  */ 

DO  10  I  =  1,  2 

CALL  ZTABLE(ALPHA(I),  ZALPHA(I)) 

10  CONTINUE 


CASE  =  0 

/*  II  IS  INDEX  FOR  N  */ 

DO  150  II  =  1,  2 

C  /*  COMPUTE  THE  DEGREE  OF  FREEDOM  FOR  CHI-SQUARE  */ 

DF  =  0. 

DO  20  I  =  1,  N( II) 

DF  =  DF  +  SAMPLE(I) 

20  CONTINUE 

DF  =  2  *  DF 


/*  FINDING  COMPONENT  TYPE  THAT  HAS  THE  MINIMUM  NUMBER  */ 
/*  OF  SAMPLE  SIZE  */ 

NSTAR  =999 
DO  25  I  =  1,  NC II) 

IF  (SAMPLE(I)  .LE.  NSTAR)  THEN 
NSTAR  =  SAMPLE(I) 

ENDIF 

25  CONTINUE 
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/*  JJ  is  INDEX  FOR  RSYS 
DO  140  JJ  =  1,  2 


*/ 


/*  COMPUTE  LAMBDA  FROM  THE  GIVEN  EQUATION  AND  */ 

/*  FIND  THE  BIGGEST  LAMBDA  */ 

BIGLAM  =  0 

DO  30  K  =  1,  N( II) 

LAMBDA(K)  =  (  -ALOG(RSYS( JJ))  /  NC II) )  /  TIME(K) 

IF(  LAMBDA(K)  .  GE.  BIGLAM)  BIGLAM  =  LAMBDA(K) 

R(K)  =  EXP(  -  LAMBDA(K)  *  TIME(K)  ) 

CONTINUE 


/*  COMPUTE  CHI-SQUARE  VALUE  FOR  GIVEN  ALPHA  AND  DF  */ 


/*  LL  IS  INDEX  FOR  ALPHA  */ 

DO  130  LL  =  1,  2 

IF  (  DF  .EQ.  1  )  THEN 
P  =  ALPHA(LL)  /  2 
CALL  ZTABLE(P,  ZALPHA(LL)) 

CHISQR  =  Z ALPHA (LL)  **  2 

ELSE  IF  (  DF  .EQ.  2  )  THEN 

CHISQR  =  -2  *  ALOG(ALPHA(LL)) 

ELSE  iF  (  DF  .GE.  3  )  THEN 

CHISQR  =  DF  *  (1  -  2. / ( 9  *  DF)  + 

ZALPHA(LL)  *  SQRT(  2. /( 9  *  DF)))  **  3 

ENDIF 


CASE  -•  CASE  +  1 
DIFF(CASE)  =  0 

DO  120  L  =  1,  NN 

LAMMAX  =  -99. 

DO  50  1=1,  N( II) 


/*  GENERATE  EXPONENTIAL  RANDOM  */ 

/*  NUMBERS  WITH  MU  =  1  */ 

CALL  LEXPN(ISEED,  B,  SAMPLE(I),  1,  0) 

T(I)  =  0 

DO  40  J  =  1,  SAMPLE(I) 

/*  CONVERT  TO  EXPONENTIAL  RANDOM  */ 
/*  NUMBERS  WITH  MU  =  LAMBDA  AND  */ 


/*  ADD  THOSE  FOR  EACH  COMPONENT  */ 
/*  TYPE  */ 


B( J)  =  B(J)  /  LAMBDA(I) 
T( I)  =  T(I)  +  B(J) 
CTIME(I.J)  =  B(J) 
CONTINUE 

LAMHAT(I)  =  SAMPLE(I)  /  T(I) 


/*  FINDING  MAXIMUM  LAMBDA  HAT  AND  ITS  */ 
/*  INDEX  */ 

IF  (  LAMHAT(I)  . GE.  LAMMAX  )  THEN 
M  =  I 

LAMMAX  =  LAMHAT(I) 

END  IF 
CONTINUE 


/*  RATIO 1  IS  FOR  WITHOUT  JACKNIFE  */ 

/*  RATI02  IS  FOR  WITH  JACKNIFE  */ 

DO  60  I  =  1,  N( II) 

RATIOl(I)  =  LAMBDA(I)  /  BIGLAM 
RATI02(I)  =  LAMHAT(I)  /  LAMMAX 
CONTINUE 


/*  PART  OF  JACKNIFE  METHOD  FOR  LAMBDA  */ 

/*  WITH  OMIT  1  VARIABLE  EACH  TIME  */ 

DO  90  I  =  1,  NC II) 

DO  80  J  =  1,  SAMPLE(I) 


LOMITf I , J)  =  0. 

DO  70  K  =  1,  SAMPLE(I) 

IF  (J  .NE.  K)  THEN 

LOMIT(I,J)  =  LOMIT(I,J)  +  CTIME(I,K) 
ENDIF 

CONTINUE 

LOMIT( I , J)  =  (SAMPLE(I)-l)  /  LOMIT(I,J) 

CONTINUE 

CONTINUE 
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/*  ADAPT  ABOVE  RESULT  TO  OUR  EQUATION  TO  */ 
/*  GET  THE  RELIABILITY  AND  TRUE  CONFIDENCE  */ 
/*  LIMIT  */ 

SUM1  =  0. 

SUM2  =  0. 

SUM3  =  0. 

SUM4  =  0. 

DO  110  I  =  1,  N('II) 


/*  NON  JACKNIFING  (  ORIGINAL  )  */ 

SUM1  =  SUM1  +  RATIOl(I)  *  T(I) 

SUM 2  =  SUM2  +  RATIOl(I)  *  TIME(I) 


/*  WITH  JACKNIFING  */ 

ROMSUM  =  0. 

DO  100  J  =  1,  NSTAR 

ROMIT(I,J)  =  LOMIT(I,J)  /  LOMIT(M,J) 
ROMSUM  =  ROMSUM  +  ROMIT(I,J) 

100  CONTINUE 

RSTAR(I)  =  NSTAR  *  RATI02(I)  - 

*  (NSTAR  -  1)  *  ROMSUM  /  NSTAR 

SUM3  =  SUM3  +  RSTAR(I)  *  T(I) 

SUM4  =  SUM4  +  RSTAR(I)  *  TIHE(I) 

110  CONTINUE 


/*  R1  ;  RELIABILITY  OF  A  SYSTEM  WITHOUT  */ 
/*  JACKNIFE  */ 
/*  R2  ;  RELIABILITY  OF  A  SYSTEM  WITH  */ 
/*  JACKNIFE  */ 


LAMHM  «  CHISQR  /  (2  *  SUM1) 

R1(L)  =  EXP(  -  LAMHM  *  SUM2  ) 

LAM1IST  =  CHISQR  /  (2  *  SUM3) 
R2(L)  =  EXP(  -  LAMHST  *  SUM4  ) 


IF  (  ABS(R1(L)  -  R2(L))  .  GE.  DIFF(CASE)  )  THEN 
DIFF(CASE)  =  ABSCRl(L)  -  R2(L)) 
DIFFRl(CASE)  =  R1(L) 

DIFFR2(CASE)  =  R2(L) 

ENDIF 

120  CONTINUE 

C  /*  NONIMSL  LIBRARY  'SHSORT'  WILL  SORT  Rl,  R2  */ 

C  /*  BY  SHELL  SORT  ALGORITHM  */ 
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CALL  SHS0RT(R1,  KEY1,  NN) 
CALL  SHS0RT(R2,  KEY2 ,  NN) 

MM  =  NN  *  (1  -  ALPHA(LL)) 

RVALl(CASE)  =  Rl(MM) 
RVAL2(CASE)  =  R2(MM) 


/*  SUBROUTINE  FINDJ  FINDS  THE  INDEX  OF  Rl,  R2  */ 
/*  WHICH  THE  VALUE  OF  IT  IS  CLOSEST  TO  RSYS.  */ 

CALL  FINDJ(R1,  NN,  RSYS(JJ),  Jl) 

CALL  FINDJ(R2,  NN,  RSYS(JJ),  J2) 

CLRl(CASE)  =  Jl  /  FLOAT(NN) 

CLR2(CASE)  =  J2  /  FLOAT(NN) 

130  CONTINUE 

140  CONTINUE 

150  CONTINUE 


WRITE(6,600) 

WRITE(6,650) 

WRITE(6,670) 

CASE  =  1 

DO  210  II  =  1,  2 
DO  200  JJ  =  1,  2 
DO  190  LL  =  1,  2 

WRITE( 6 ,700)  N( II),  RSYS(JJ),  ALPHA (LL) ,  RVALl(CASE), 
*  RVAL2(CASE) ,  DIFF(CASE) ,  CLRl(CASE),  CLR2(CASE) 

WRITE(6,888)  DIFFRl(CASE) ,  DIFFR2(CASE) 

CASE  =  CASE  +  1 

190  CONTINUE 

200  CONTINUE 

WRITE( 6 , 777)  (  SAMPLE(J),  J=1,N(II)  ) 

WRITE(6,999)  (  R(J),  J=1,I  ,11)  ) 

210  CONTINUE 


600  FORMAT('l' ,5(/).7X,'****  RELIABILITY  OF  SERIES  SYSTEM  ****’) 

650  F0RMAT(///,T50, 'Rl,  CL1  ;  WITHOUT  JACKNIFING', 

*  /,T50.'R2,  CL2  ;  WITH  JACKNIFING') 

670  F0RMAT(///,T6, 'NUMBER  OF' ,T19 , 'RELIABILITY' ,T33, 'ALPHA' ,T46,'Rl'  , 

*  T56. ’ R2 ' ,T63, 'MAX  (Rl  -  R2) ' ,T79, 'TRUE' ,T87  'TRUE' ,/, 

*  T6,' COMPONENTS ' ,T1 9, 'OF  SYSTEM’ ,T79, 'C. L.  lf,T87,'C.L.  2*, 

*  / ,T5 ,89( ' -' )) 
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700  F0RMAT(/,T8,I5,T22,F5.  3,T33,F5.  2,T43,F7. 4,T53,F7. 4,T65,F8. 4,T78, 
*  F7. 4,T86,F7.  4) 


777  F0RMAT(/,T3.' SAMPLE  SIZE  ARE  5(2X,I3)  ) 
888  FORMAT(T60,  Rl=' ,F6. 4, ’  R2=',F6.4) 

999  F0RMAT(/,T3, 'R(I)  ;  ' ,5(1X,F5. 3)) 

STOP 

END 


*  * 
ititititititititMtititititititititMtitititititititititititMtMtitititititMtitititititititMtMcititititititMtititicitititititititit 


SUBROUTINE  ZTABLE ( ALPHA, ZALPHA) 

«  SUBROUTINE  ZTABL2  COMPUTES  RIGHT  PERCENT  POINT  ZALPHA  FROM  » 
«  RIGHT  CUMULATIVE  PROBABILITY  ALPHA  » 

REAL  ALPHA,  ZALPHA 

IF  ((ALPHA  .GT.  0.0)  .OR.  (ALPHA  .  LT.  1))  THEN 

W  =  -  AL0G(4  *  ALPHA  *  (  1  -  ALPHA  )) 

ZALPHA  =  SQRT(W  *  (2.06118  -  (  5.72622  /  (W  +  11.6406)))) 

IF  (ALPHA  .GT.  0.5)  THEN 
ZALPHA  =  -  ZALPHA 

END  IF 

ENDIF 

RETURN 

END 


************************************************************************ 
it  it 

ititititititititititititititititititititititiiitititititititititititiiiii.ititititititititititititititititititititititititititititititititititit 


SUBROUTINE  SHSORT(A,KEY,N) 
DIMENSION  A(N) ,KEY(N) 

Ml=l 

6  M1=M1*2 

IF  (Ml  .LE.  N)  GO  TO  6 
Ml=Ml/2-l 
MM=MAX0(Ml/2, 1) 

GO  TO  21 

20  MM=MM/2 

IF  (MM  .LE.  0)  GO  TO  100 

21  K=N-MM 
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22  DO  1  J=1,K 
II=J 

11  IM=II+MM 

IF  (A( IM)  .GE.  A(II))  GO  TO  1 
TEMP=A(II) 

IT=KEy(II) 

A( II)=A( IM) 

KEY( II)-KEY( IM) 

A(IM)=TEMP 
KEY( IM)=IT 
II=II-MM 

IF  (II  .GT.  0)  GO  TO  11 
1  CONTINUE 
GO  TO  20 
100  RETURN 
END 


*  * 

SUBROUTINE  FINDJ(A,  NN,  R,  J) 

«  SUBROUTINE  FINDJ  FINDS  THE  INDEX  OF  ARRAY  A  WHICH  THE  » 
«  VALUE  OF  IT  IS  CLOSEST  TO  R.  » 

REAL  A(NN) ,  R,  VALUE 
INTEGER  J 

VALUE  =  ABS(A(NN)  -  R) 

DO  100  I  =  NN" 1 j  1,  -1 

IF  (ABS(A( I)  -  R)  .LT.  VALUE)  THEN 
VALUE  =  ABS(A(I)  -  R) 

ELSE 

J  =  I  +  1 
RETURN 

ENDIF 

100  CONTINUE 

RETURN 
END 
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APPENDIX  C.  FORTRAN  CODE  FOR  INTERVAL  ESTIMATION 
PROCEDURE  -  WEIBULL  CASE 


c 

c 

c 

c 

c 

c 

c 


PROGRAM  WEIBUL 

******5V****5V***********rt>V*rt*******rt*rt*rt5V***foVjV*A*A*****rt***&?V***** 
it  it 

*  THIS  PROGRAM  COMPUTES  THE  RELIABILITY  OF  A  SERIES  SYSTEM  * 

*  USING  TWO  DIFFERENT  METHODS  FOR  THE  CASE  WHEN  FAILURE  * 

*  OF  EACH  COMPONENTS  IS  WE I BULLY  DISTRIBUTED.  * 

*  AND  ALSO  COMPARES  THE  ESTIMATE  OF  SHAPE  PARAMETERS  (  BETA  )  * 


c 

* 

WHICH  IS  COMPUTED  FROM  TWO  DIFFERENT  METHODS. 

* 

c 

it 

* 

c 

* 

THE  SERIES 

SYSTEM  CONSIDERED  IN  THIS  PORGRAM  HAS 

* 

c 

it 

N  COMPONENT  TYPES  AND  EACH  COMPONENT  HAS  SAME  NUMBER  OF 

* 

c 

p 

* 

"if 

SAMPLE  SIZES. 

* 

it 

c 

p 

•ff 

BELOWS  ARE 

GIVEN  OR  ASSUMED  ; 

it 

it 

c 

* 

BETA 

SHAPE  PARAMETER  OF  WEIBULL  DISTRIBUTION 

it 

c 

it 

ALPHA 

SIGNIFICANCE  LEVEL 

* 

c 

it 

X 

TEST  TIME 

* 

c 

P 

it 

RSYS 

RELIABILITY  OF  A  SERIES  SYSTEM 

* 

it 

c 

p 

it 

THESE  ARE  VARIABLES  USED  ; 

it 

it 

c 

it 

BETHAT 

BETA  HAT  FROM  THE  FORMULA 

it 

c 

iV 

BETNEW 

BETA  HAT  FROM  NEWTON  -  RAPHSON  METHOD 

it 

c 

* 

CHISQR 

CHI-SQUARE  VALUE  FOR  GIVEN  ALPHA  AND  DF 

it 

c 

Vr 

CLR 

TRUE  CONFIDENCE  LIMIT  FROM  THE  FORMULA 

it 

c 

* 

CLRNEW 

TRUE  CONFIDENCE  LIMIT  FROM  THE  N-R  METHOD 

it 

c 

EXPO 

EXPONENTIAL  RANDOM  NUMBERS  WITH  MU=1 

it 

c 

* 

KSUBM 

K  SUB  M  VALUE 

* 

c 

* 

LAMBDA 

SCALE  PARAMETER  OF  WEIBULL  DISTRIBUTION 

it 

c 

LAMBIG 

MAXIMUM  VALUE  OF  LAMBDA  FROM  THE  N-R  METHOD 

it 

c 

* 

LAMHAT 

ESTIMATES  OF  LAMBDA  FROM  THE  FORMULA 

it 

c 

* 

LAMMAX 

MAXIMUM  VALUE  OF  LAMBDA  FROM  THE  FORMULA 

it 

c 

* 

LAMMU 

FINALLY  COMPUTED  LAMBDA  FROM  THE  FORMULA 

it 

c 

* 

LAMSTR 

LAMBDA  TO  THE  BETA 

it 

c 

JU 

LAMNEW 

ESTIMATES  OF  LAMBDA  FROM  THE  N-R  METHOD 

it 

c 

* 

LAMUNR 

FINALLY  COMPUTED  LAMBDA  FROM  THE  N-R  METHOD 

it 

c 

* 

RATIO 

RATIO  OF  LAMBDA  FOR  THE  FORMULA 

it 

c 

?v 

RATNEW 

RATIO  OF  LAMBDA  FOR  THE  N-R  METHOD 

it 

c 

* 

RHAT 

R(  500  *  (1-ALPHA)  )  FOR  THE  FORMULA 

it 

c 

Vr 

RHTNEW 

R(  500  *  (1-ALPHA)  )  FOR  THE  N-R  METHOD 

it 

c 

XBENEW 

X  TO  THE  BETNEW 

it 

c 

* 

XHAT 

SUM  OF  XIJHAT 

it 

c 

* 

XIJHAT 

W  TO  THE  BETHAT 

it 

c 

* 

XIJNEW 

W  TO  THE  BETNEW 

it 

c 

* 

XNEW 

SUM  OF  XIJNEW 

it 

c 

* 

XTOBET 

X  TO  THE  BETHAT 

it 

c 

W.WEIB 

WEIBULL  RANDOM  NUMBERS 

it 
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*  ZALPHA  ;  RIGHT  PERCENTILE  POINT  (NORMAL  DISTRIBUTION)  * 

*  iV 
*rt***-.Wrrt***rt**iV*VoW«VA**&iV******A***5WiiV!ViY*?VAAiV***********iV**ft****y«V 


PARAMETER  (  N=5,  M=30  ,  NN=  500  ) 

PARAMETER  (  N=15,  M=30  ,  NN=  500) 

REAL  BETA(N) ,  LAMBDA(N) ,  X(N)S  XHAT(N) ,  XTOBET(N),  XX(N,M) 
REAL  W(N,M) ,  Y(N,M) ,  BETHAT(N) ,  WBHAT(N,M) ,  LAMHAT(N) 

REAL  RATIO(N) ,  LAMSTR(N) ,  LAMMAX,  RSYS(2),  KSUBM(N) ,  LAMMU 
REAL  EXPO(M) ,  ALPHA(2) ,  ZALPHA(2),  RHAT(10),  RSL(NN) ,  CHISQR 
REAL  B(M\  KEY(M) ,  KEYl(NN) ,  KEY2(NN) ,  XIJHAT(N,M) 

REAL  WEIB(M),  BETNEW(N) ,  XBENEW(N),  XNEW(N),  XIJNEW(N,M) 

REAL  LAMNEW(N) ,  LAMBIG,  RATNEW(N) ,  LAMUNR,  RSLNEW(NN) 

REAL  RHTNEW(IO),  CLR(IO),  CLRNEW(IO) 

REAL  BHTBAR ,  BNRBAR,  BHTMSE,  BNRMSE 

INTEGER  S(N) ,  DF,  ICOUNT,  SAMPLE(N),  CASE 

DATA  ISEED  /  1736  / 

DATA  ALPHA  /  .2,  .05  / 

DATA  RSYS  /  .  90,  .  975  / 

DATA  BETA  /  1.  2,1.  4,1.  6,1.  8,1.  5  / 

DATA  BETA  /l.  2 , 1.  4, 1.  6 , 1.  8, 1.  5 , 1.  3, 1.  7 , 1.  5 , 1.  2, 1.  6, 1.  9 , 1.  8, 1.  4, 

*  1.6,1.  1/ 

DATA  BETA  /  2.  2,2. 4,2. 6,2. 8,2. 5  / 

DATA  BETA  /2.  2,2.  4,2.  6,2.  8,2.  5,2.  3,2.  7,2.  5,2.  2,2.  6,2.  9,2.  8,2.  4, 

*  2.6,2.  1/ 

DATA  X  /  2,3,7,8,10/ 

DATA  X  /2, 3, 7, 8, 10, 3, 7, 10, 1,7, 8, 3, 10, 1,8/ 

DATA  SAMPLE  /9, 7, 10, 8, 6/ 

DATA  SAMPLE  /30, 21, 10, 15,26/ 

DATA  SAMPLE  /  30,63,75,98,26  / 

DATA  SAMPLE  /6, 7 ,5 ,6 , 7 ,8 ,9 ,5 ,8,7 ,9 , 10, 7 ,9 ,6/ 

DATA  SAMPLE  /3, 7, 10, 15,20, 15 ,7,5,30,20, 10, 7, 13,21,30/ 

DATA  SAMPLE  /15, 40, 35, 17, 26, 67, 50, 65, 80, 32, 95, 100, 15, 45, 30/ 

DATA  SAMPLE  /  15  *  15  / 

CASE  =  1 

DO  10  I  =  1,  N 

IF  (  SAMPLE(I)  .LE.  15)  KSUBM(I)  =  1.40 
IF  (  SAMPLE(I)  .LE.  30)  KSUBM(I)  =  1.50 
10  CONTINUE 


DO  220  JJ  =  1,  2 

C  /*  FINDING  THE  LAMBDA  AND  LAMBDA  STAR  FROM  THE  GIVEN  DATA  */ 

DO  20  I  =  1,  N 
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LAMBDA(I)  =  ((  -  ALOG(RSYS(JJ))  /  N)  **  (1. /BETA(I)))/  X(I) 
LAMSTR(I)  =  LAMBDA(I)  **  BETA(I) 

S(I)  =0.84  *  SAMPLE(I)  -  0.5 
20  CONTINUE 

C  /*  SUBROUTINE  ZTABLE  COMPUTES  THE  RIGHT  PERCENT  POINT  ZALPHA  */ 

C  /*  FROM  RIGHT  CUMULATIVE  PROBABILITY  ALPHA  */ 

DO  40  I  =  1,  2 

CALL  ZTABLE (ALPHA( I),  ZALPHA(I)) 

40  CONTINUE 


C  /*  COMPUTE  THE  DEGREE  OF  FREEDOM  FOR  CHI-SQUARE  */ 

DF  =  0 

DO  60  I  =  1,  N 

DF  =  DF  +  SAMPLE(I) 

60  CONTINUE 

DF  =  2  *  DF 


C  /*  COMPUTE  CHI-SQUARE  VALUE  FOR  GIVEN  ALPHA  AND  DF  */ 

DO  200  LL  =  1,  2 

IF  (  DF  .EQ.  1  )  THEN 
P  =  ALPHA (LL)  /  2 
CALL  ZTABLECP,  ZALPHA(LL)) 

CHISQR  =  ZALPHA(LL)  **  2 

ELSE  IF  (  DF  .EQ.  2  )  THEN 

CHISQR  =  -2  *  ALOG(ALPHA(LL)) 

ELSE  IF  (  DF  .GE.  3  )  THEN 

CHISQR  =  DF  *  (1  -  2./(9  *  DF)  + 

*  ZALPHA(LL)  *  SQRT(  2./(9  *  DF)))  **  3 

ENDIF 

LAMMAX  =  0 

LAMB  I G  =  0 

DO  180  ITER  =  1,  NN 
DO  140  I  =  1,  N 

/*  GENERATE  EXPONENTIAL  RANDOM  NUMBER  WITH  MU=1  */ 
CALL  LEXPNC I SEED,  EXPO,  SAMPLE(I),  1,  0) 

/*  CONVERT  EXPONENTIAL  RANDOM  NUMBER  WITH  MU=1  */ 
/*  TO  MU=LAMBDA  AND  GET  WEIBULL  RANDOM  NUMBERS  */ 
DO  80  J  =  1,  SAMPLE(I) 

XX(I,J)  =  EXPO(J)  /  LAMSTR(I) 

W(I,J)  =  XX(I,J)  **  (l/BETACI)) 

Y( I , J)  =  ALOG(W(I, J)) 

WEIB(J)  =  W(I,J) 
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B(J)  =  Y(I,J) 

CONTINUE 

/*  SUBROUTINE  SHSORT  WILL  SORT  ARRAY  B  IN  */ 

/*  ASCENDING  ORDER  BY  SHELL  SORT  ALGORITHM  */ 

CALL  SHSORT( B ,  KEY,  SAMPLE(I)) 

/*  FINDING  BETA  HAT  (BETHAT)  BY  THE  FORMULA  */ 

SUM1  =  0 
SUM2  =  0 

DO  100  J  =  1,  S(I) 

SUM1  =  SUM1  +  B(J) 

CONTINUE 

DO  120  J  =  S(I)+1,  SAMPLE(I) 

SUM2  =  SUM2  +  B(J) 

CONTINUE 

BETHAT(I)  =  (SAMPLE(I))  *  KSUBM(I)  / 

(S(I)*SUM2  /  (SAMPLE(I)  -  S(I))  -  SUM1) 

TEMP  =  BETA(I) 

/*  SUBROUTINE  NEWTON  WILL  COMPUTE  THE  ESTIMATE  */ 
/*  OF  BETA  (BETNEW)  BY  NEWTON -RAPHSON  METHOD  */ 
CALL  NEWTON(WEIB,  SAMPLE(I),  TEMP,  BETNEW(I)) 

XTOBET(I)  =  X(I)  **  BETHAT(I) 

XBENEW(I)  =  X(I)  **  BETNEW(I) 

XHAT(I)  =  0. 

XNEW(I)  =  0. 

DO  130  J  =  1,  SAMPLE(I) 

XIJHAT(I,J)  =  W(I,J)  **  BETHAT(I) 

XHAT(I)  =  XHAT(I)  +  XIJHAT(I,J) 

XIJNEW( I , J)  =  W(I,J)  **  BETNEW(I) 

XNEW(I)  =  XNEW(I)  +  XIJNEWC I , J) 

CONTINUE 

LAMHAT(I)  =  SAMPLE(I)  /  XHAT(I) 

LAMNEW(I)  =  SAMPLE(I)  /  XNEW(I) 

IF  (LAMHAT(I)  . GT.  LAMMAX  )  THEN 
LAMMAX  =  LAMHAT(I) 

END  IF 

IF  (LAMNEW(I)  .GT.  LAMBIG  )  THEN 
LAMBIG  =  LAMNEW(I) 

ENDIF 

CONTINUE 


SUM3N  =  0 
SUH4  =  0 
SUM4N  =  0 


DO  160  I  =  1,  N 

RATIO(I)  =  LAMHAT(I)  /  LAMMAX 
RATNEW(I)  =  LAMNEW(I)  /  LAMBIG 
SUM3  =  SUM3  +  RATIO(I)  *  XHAT(I) 

SUM4  =  SUM4  +  RATIO(I)  *  XTOBET(I) 

SUM3N  =  SUM3N  +  RATNEW(I)  *  XNEW(I) 
SUM4N  =  SUM4N  +  RATNEW(I)  *  XBENEW(I) 

160  CONTINUE 

LAMMU  =  CHISQR  /  (  2  *  SUM3  ) 

RSL(ITER)  =  EXP(  -  LAMMU  *  SUM4  ) 

LAMUNR  =  CHISQR  /(  2  *  SUM3N  ) 
RSLNEW(ITER)  =  EXP(  -  LAMUNR  *  SUM4N  ) 

180  CONTINUE 

CALL  SHSORT(RSL,  KEY1,  NN) 

CALL  SHSORT(RSLNEW,  KEY2,  NN) 

KK  =  (  1  -  ALPHA(LL))  *  NN 


/*  RHAT  IS  THE  RELIABILITY  OF  THE  SERIES  SYSTEM  COMPUTED  */ 
/*  BY  THE  FORMULA  */ 

/*  RHTNEW  IS  THE  RELIABILITY  OF  THE  SERIES  SYSTEM  COMPUTED  */ 
/*  BY  THE  NEWTON  -  RAPHSON  METHOD  */ 

RHAT( CASE)  =  RSL(KK) 

RHTNEW(CASE)  =  RSLNEW(KK) 

/*  SUBROUTINE  FINDJ  WILL  FIND  THE  TRUE  CONFIDENCE  LIMIT  */ 
CALL  FINDJ(RSL,  NN,  RSYS(JJ),  Jl) 

CALL  FINDJ(RSLNEW,  NN,  RSYS(JJ),  J2) 


CLR(CASE)  =  Jl  /  FLOAT(NN) 

CLRNEW(CASE)  =  J2  /  FLOAT(NN) 

CASE  =  CASE  +  1 

200  CONTINUE 


220  CONTINUE 

DO  240  I  =  1,  N 

WRITE(6,500)  BETA(I),  SAMPLE(I) 
240  CONTINUE 

WRITE(6,550)  N,  NN 
WRITE(6,600) 
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WRITE(6,650) 

WRITE(6,670) 

CASE  =  1 

DO  280  JJ  =  1,  2 

DO  260  LL  =  1,  2 

WRITE ( 6 , 7 00 )  N , RSYS ( JJ) , ALPHA ( LL) , RHAT( CASE ) , RHTNEW( CASE ) , 
*  CLR(CASE) ,  CLRNEW(CASE) 

CASE  =  CASE  +  1 

260  CONTINUE 
280  CONTINUE 


500  FORMAT(T10,' INITIAL  BETA  : ' ,3X,F4. 1, 

*  T45, 'NUMBER  OF  SAMPLE  SIZES  :',2X,I3) 

550  F0RMAT(T10, 'NUMBER  OF  COMPONENTS  : ' , 2X, 13 , 

*  /,T10, 'NUMBER  OF  REPLICATIONS  :',2X,I3) 

600  FORMATCl'  ,5(/),7X,'****  RELIABILITY  OF  SERIES  SYSTEM  ****') 

650  FORMAT( /// ,T35 , 'Rl,  CL1  ;  USING  THE  FORMULAR* , 

*  /,T35 ,'R2,  CL2  ;  USING  THE  NEWTON  RAHPSON  METHOD') 

670  F0RMAT(///,T6,  NUMBER  OF' .T19, 'RELIABILITY' ,T33, 'ALPHA' ,T46,'Rl' , 

*  T56 ,'R2' ,T66, 'TRUE  ,T74,'TRUE' ./, 

*  T6, Components' ,ti9, 'of  sYSTEMf ,T66, 'c. l.  i',T74,'c.l.  2', 

*  / >T5 , 76( ' - ' ) ) 

700  FORMAT( / ,T8, 15 ,T22 ,F5.  3^33^5.2^43^7.4^53^7.4^65, 

*  F7. 4,T73 ,F7.  4) 


STOP 

END 


****■ 

* 

***•>’;■ 


'*******rt**rty?*******rt***rt***Art****rt-.V*yfrt**rt*****rtArtiWrrt***rt*** 

* 


•iticicMcMtMt  A*A**A**AAAA*A**yfA*A*iVAA*A**5Wf*A***A 


SUBROUTINE  ZTABLE( ALPHA, ZALPHA) 

«  SUBROUTINE  ZTABL2  COMPUTES  RIGHT  PERCENT  POINT  ZALPHA  FROM  » 
«  RIGHT  CUMULATIVE  PROBABILITY  ALPHA  » 

REAL  ALPHA,  ZALPHA 


IF  ((ALPHA  .GT.  0.0)  .OR.  (ALPHA  .  LT.  1))  THEN 
W  =  -  AL0G(4  *  ALPHA  *  (  1  -  ALPHA  )) 

ZALPHA  =  SQRT(W  *  (2.06118  -  (  5.72622  /  (W  +  11.6406)))) 

IF  (ALPHA  .  GT.  0.5)  THEN 
ZALPHA  =  -  ZALPHA 

END  IF 
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END  IF 


RETURN 

END 


iWnV*)V*VoViWf>V*iV5Wr**iV5WfV()V5V*)!r;V*5Wf*5WfiV)V>V',Wr5'nWr»VVriV*5ViV)V5V**A'5VA*y('A*yryr5V*A*>VAA'5V5'r?Viln^* 
*  * 
iV*5V**5V**Vf**iV*?V*A***5VA***5V***')V>ViV>V?V****5V*5ViV*Vr5V**iV***>V*A>V*>V>V*5V*»VA**>V*>'nViV*** 


SUBROUTINE  SHSORT(A,KEY,N) 
DIMENSION  A(N) ,KEY(N) 

Ml=l 

6  M1=M1*2 

IF  (Ml  .IE.  N)  GO  TO  6 
Ml=Ml/2-l 
MM=MAX0(Ml/2, 1) 

GO  TO  21 

20  MM=MM/2 

IF  (MM  .LE.  0)  GO  TO  100 

21  K=N-MM 

22  DO  1  J=1,K 
II*J 

11  IM-II+MM 

IF  (A(IM)  .GE.  A(II))  GO  TO  1 
TEMP=A( II) 

IT=KEY( II) 

A(II)=A(IM) 

KEY(II)=KEY(IM) 

A( IM)=TEMP 
KEY(IM)=IT 
II=II-MM 

IF  (II  .GT.  0)  GO  TO  11 
1  CONTINUE 
GO  TO  20 
100  RETURN 
END 


Vf***Vf*Vf**************yf******Vf**Vf***A****vr*Vfyfyf*************Vf**yf**yf****yf** 
*  * 


SUBROUTINE  NEWTON(TIME,  M,  BETA,  BETNEW  ) 

«  SUBROUTINE  NEWTON  WILL  COMPUTE  THE  ESTIMATE  OF  BETA  (BETHAT)  » 
«  BY  NEWTON  -  RAPHSON  METHOD  » 

REAL  TIME(M) ,  LNTIME(IOO),  A(3),  LNTS 
INTEGER  R 

C  DATA  TS  /  6  / 

C  =0 
ITER  =  0 
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R  =  M 
LNTS  =  ALOG(TS) 


DO  20  I  =  1,  R 

LNTIME(I)  =  ALOG(TIME(I)) 
C  =  C  +  LNTIME(I) 

20  CONTINUE 


C  =  C  /  R 

30  DO  60  J  =  1,  3 
SUM  =  0 

DO  40  K  =  1,  R 

IF  (  J  .EQ.  1  )  THEN 

SUM  =  SUM  +  TIME(K)  **  BETA 

ELSE 

SUM  =  SUM  +  ((TIME(K)  **  BETA ) * ( LNTI ME ( K )  **  (J-l))) 
ENDIF 

40  CONTINUE 

AC J)  =  SUM 
60  CONTINUE 

/*  FUNCTION  FPRIME  IS  THE  DERIVATIVES  OF  FUNCTION  F  */ 

QUOT  =  A(2)  /  A( 1) 

FPRIME  =  A(3)  /  A(  1)  -  (QU0T**2)  +  ( ( 1/BETA)**2) 

F  =  QUOT  -  (1  /  BETA)  -  C 

/*  BETA  IS  UPDATED  EACH  TIME  AND  CHECK  IF  IT  CONVERGES  */ 

BETA  =  BETA  -  F  /  FPRIME 

ITER  =  ITER  +  1 

IF  (BETA  ,GT.  25.)  GOTO  100 

IF  (  ABS(F)  .GT.  0.0001  )  GOTO  30 

ALPHA  =  (A( 1)  /  R)  **  (1/BETA) 

BETNEW  =  BETA 


RETURN 

100  WRITE(6,*)  'DID  NOT  CONVERGE' 

WRITE(6,*)  'TRY  AGAIN  WITH  BETTER  ESTIMATE  OF  BETA' 

RETURN 

END 


*********************************************************?Wf*****Vf***** 

* 

******?Wt****A**Vf*****5Hr******************A*Vf  ************************** 
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SUBROUTINE  FINDJ(A,  NN,  R,  J) 

«  SUBROUTINE  FINDJ  FINDS  THE  INDEX  OF  ARRAY  A  WHICH  THE  » 
«  VALUE  OF  IT  IS  CLOSEST  TO  R.  » 

REAL  A(NN) ,  R,  VALUE 
INTEGER  J 

VALUE  =  ABS(A(NN)  -  R) 

DO  100  I  =  NN-1,  1,  -1 

IF  (ABS(A( I)  -  R)  .LT.  VALUE)  THEN 
VALUE  =  ABS(A( I)  -  R) 

ELSE 

J  =  I  +  1 
RETURN 

ENDIF 

100  CONTINUE 

RETURN 

END 


o  o  o  o  ooooooooooooo 


APPENDIX  D.  FORTRAN  CODE  FOR  INTERVAL  ESTIMATION 
PROCEDURE  -  NORMAL  CASE 

PROGRAM  NORMAL 


******A**A****AA******A*iV*rt**>V5V***W*rt>V*yf*****>Y**>V**>V****>Wf*iVft>V**ftiV 


*  * 

*  THIS  PROGRAM  DETERMINE  THE  ACCURACY  OF  AN  APPROXIMATE  * 

*  LOWER  CONFIDENCE  BOUND  FOR  P(  X  >  Y  ).  * 

*  * 

*  PROGRAM  NORMAL  IS  THE  CASE  WHEN  Y  IS  GIVEN  A  VALUE  YO  * 

*  &  X  IS  NORMALLY  DISTRIBUTED  WITH  * 

*  UNKNOWN  PARAMETERS.  * 

*  * 


*****rtiY********rt*******A*iYiYiYiY**ynYiVArtiWrA****iY-.Y*rt*iYrtrt**yr***yn'n'nYyoY*ynY 


REAL  YO,  P,  ZP,  MUX,  SUMX,  SUMX2 ,  XBAR 

REAL  TEMPI,  TEMP2,  TEMP3,  SIGHAT 

REAL  R(3),  SIGMAX(3) ,  ALPHA(2) ,  X(100),  Y(100) 

REAL  RLl(54,i000) ,  ARLl(lOOO),  BRL1(54),  KEY(IOOO),  CIRL1(54) 
REAL  T1RL1(54),  T1CI1(54),  T2RL1(54),  T2CI1(54) 

REAL  TRL1( 2) ,  TCI 1(2) 

REAL  ZHAT(IOOO),  ZKNIFE(S4),  ZVAR(54),  ZBAR(54) 

REAL  MUJUNK(54) 

INTEGER  NUMX( 3) ,  NUMY(3),  CASE,  LINE 


DATA  R  /. 95,. 99,. 995/,  SIGMAX  /. 5 , 1 ,20/ ,  NUMX  /10,25,75/ 

DATA  NUMY  /10,25,75/  ,  ALPHA  /  .2,  .1/ 

DATA  ISEED,  JSEED  /  4875,7981  / 

DATA  YO  /  400  / 

CASE  =  0 

/*  II  IS  THE  INDEX  FOR  R.  R(l)  =  .95,  R(2)  =  .99,  R(3)  =  .995  */ 

DO  500  II  =  1,  3 


/*  P  IS  THE  RIGTH( UPPER)  CUMULATIVE  PROBABILITY  */ 

P  =  1  -  R( II) 


/*  SUBROUTINE  ZTABL2  COMPUTES  RIGHT  PERCENT  POINT  ZP  */ 

/*  FROM  RIGHT  CUMULATIVE  PROBABILITY  P  */ 

CALL  ZTABL2(  P,  ZP  ) 


C  /*  JJ  IS  THE  INDEX  FOR  SIGMA  OF  X.  */ 

C  /*  SIGMAX(l)  =  .5,  SIGMAX( 2)  =  1,  SIGMAX(3)  =  3  */ 
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DO  400  JJ  =  1,  3 

MUX  =  ZP  *  SIGMAX(JJ)  +  YO 


/*  KK  IS  THE  INDEX  OF  NUMBER  OF  X.  */ 

/*  NUMX(l)  =  10,  NUMX(2)  =  25,  NUMX(3)  =  75  */ 

DO  300  KK  =  1,  3 

C  /*  LL  IS  THE  INDEX  FOR  ALPHA.  */ 

C  /*  ALPHA ( 1)  =  .2,  ALPHAC2)  =  .1  */ 

DO  200  LL  =  1,  2 


CALL  ZTABL2( ALPHA ( LL) ,  ZALPHA) 

CASE  =  CASE  +  1 

TEMP  =  (NUMX(KK)  -  1.  )  /  NUMX(KK) 
TEMPO  =  SQRT(TEMP) 

MUJUNK(CASE)  =  MUX 


C  /*  REPLICATE  1000  TIMES  FOR  EACH  CASES  */ 

DO  100  I  =  1,  1000 

C  /*  USE  NORMAL  RANDOM  NUMBER  GENERATOR  TO  GET  */ 

C  /*  NUMX(KK)  NUMBER  OF  X.  */ 

CALL  LNORM(ISEED,X,NUMX(KK) ,2,0) 

C  /*  THIS  PART  IS  TO  GET  SAMPLE  MEAN(XBAR)  AND  */ 

C  /*  SAMPLE  VARIANCE(XVAR)  */ 


DO  50  MM  =  1,  NUMX(KK) 

X(MM)  =  SIGMAX(JJ)  *  X(MM)  +  MUX 
50  CONTINUE 


C  /*  SUBROUTINE  VAR  WILL  COMPUTE  SAMPLE  MEAN  AND  */ 

C  /*  SAMPLE  VARIANCE  */ 

CALL  VAR(X,  NUMX(KK) ,  XBAR,  XVAR) 

C  /*  NOW  WE  COMPUTE  THE  LOWER  CONFIDENCE  BOUND  */ 

TEMPI  =  (XBAR  -  YO)  /  SQRT(XVAR  *  TEMPO) 

TEMP2  =  1.  /  NUMX(KK)  + 

*  ((XBAR  -  YO)  **  2)  /  ( 2*( NUMX( KK)+1 )*  XVAR) 

TEMPI  =  TEMPI  -  ZALPHA  *  SQRT(TEMP2)  *  TEMP 


C  /*  SUBROUTINE  ZTABL1  COMPUTE  RIGHT  CUMULATIVE  PROBA-  */ 

C  /*  BILITY  FROM  RIGHT  PERCENT  POINT.  */ 

CALL  ZTABL1 (TEMPI,  ARL1(I)) 
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ARL1(I)  =  1.  -  ARL1(I) 

CALL  JKNIFE(X,  NUMX(KK) ,  YO,  ZHAT(I)) 

100  CONTINUE 

CALL  VARCZHAT,  1000,  ZBAR(CASE) ,  ZVAR(CASE)) 

SIGHAT  =  1.  /NUMX(KK)  +  ZBAR(CASE)**2  /  (2*(NUMX(KK)+1)) 

SIGHAT  =  SQRT( SIGHAT)  *  TEMP 

SIGHAT  =  ZBAR(CASE)  /  TEMPO  -  SIGHAT  *  ZALPHA 

CALL  ZTABLK SIGHAT,  ZKNIFE(CASE) ) 

ZKNIFE(CASE)  =  1  -  ZKNIFE(CASE) 


C  /*  NON-IMSL  LIBRARY  'SHSORT'  WILL  SORT  ARL1  BY  SHELL  */ 

C  /*  SORT  ALGORITHM.  */ 

CALL  SHS0RT(ARL1,  KEY,  1000) 

DO  150  I  =  1,  1000 
RL1(CASE,I)  =  ARL1(I) 

150  CONTINUE 


C  /*  SUBROUTINE  FINDJ  FINDS  THE  INDEX  OF  ARL20  WHICH  THE  */ 

C  /*  VALUE  OF  IT  IS  CLOSEST  TO  R.  */ 

CALL  FINDJ(ARL1,  R( II) ,  J) 


C  /*  DIVIDING  THE  INDEX  BY  1000  WILL  GIVE  US  TRUE  */ 

C  /*  CONFIDENCE  LEVEL.  */ 

CIRLl(CASE)  =  J  /  1000. 

MM  =  1000  *  (1  -  ALPHA(LL)) 

BRLl(CASE)  =  RLl(CASE.MM) 

CALL  TNYVAL(R(II),  SIGMAX(JJ),  NUMX(KK) , 

*  ALPHA ( LL) ,  ZALPHA,  MUX,  YO,  TRL1,  TCI1) 

TlRLl(CASE)  =  TRLl(l) 

TlCIl(CASE)  =  TCIl(l) 

T2RL3 '  ’ASE)  =  TRL1(2) 

T2CI1(CASE)  =  TCI1(2) 


200  CONTINUE 

300  CONTINUE 
400  CONTINUE 
500  CONTINUE 
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WRITE(6,600) 

WRITE(6,630) 

1  =  1 
LINE  =  0 

DO  550  II  =  1,  3 
DO  550  JJ  =  1,  3 
DO  550  KK  =  1,  3 

IF  (LINE  .GE.  3)  THEN 
WRITE(6,800) 

LINE  =  0 
ENDIF 


* 

* 

* 


WRITE(6,710)  RC II) ,  SIGMAX(JJ),  NUMX(KK), 

BRL1(I) ,  BRLKI+l),  CIRLl(I) ,  CIRL1(I+1), 
ZKNIFE(I),  ZKNIFECI+1), 

^  MUJUNK(I),  MUJUNK( 1+1) 


LINE  =  LINE  +  1 
550  CONTINUE 


WRITE(6,610) 

WRITE(6,650) 

1  =  1 
LINE  =  0 

DO  560  II  =  1,  3 
DO  560  JJ  =  1,  3 
DO  560  KK  =  1,  3 

IF  (LINE  ,  GE.  3)  THEN 
WRITE(6,800) 

LINE  =  0 
ENDIF 

WRITE(6, 700)  R(II) ,  SIGMAX(JJ),  NUMX(KK) , 

j.  _  +  TlRLl(I),  T1RLKI+1),  TlCIl(I),  T1CI1(I+1) 

LINE  =  LINE  +  1 
560  CONTINUE 


WRITE(6,620) 

WRITE(6,650) 

I  =  1 
LINE  =  0 

DO  570  II  =  1,  3 
DO  570  JJ  =  1,  3 
DO  570  KK  =  1,  3 

IF  (LINE  .GE.  3)  THEN 
WRITE (6, 800) 

LINE  =  0 
ENDIF 
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WRITE(6, 700)  R( II) ,  SIGMAX(JJ),  NUMX(KK) , 

*  T2RL1(I),  T2RL1(I+1),  T2CI1(I),  T2CI1(I+1) 

1  =  1  +  2 
LINE  =  LINE  +  1 
570  CONTINUE 


600  FORMATC'l' ,5(/) ,7X, *****  Y  IS  GIVEN  A  VALUE  OF  Y0  ****') 

610  FORMATC'l '  ,5(/),7X,'****  Y  IS  GIVEN  A  VALUE  OF  Y0  ****', 

*  //.20X,  ~  X  IS  TRUNCATED  NORMAL  WITH  PROBABILITY  .95  — ') 
620  FORMATC'l' ,5(/) ,7X, '****  Y  IS  GIVEN  A  VALUE  OF  YO  ****', 

*  // ,20X, 1 --  X  IS  TRUNCATED  NORMAL  WITH  PROBABILITY  .90  — ') 

630  FORMAT(3(/),5X,'  CASE  ’ ,6X,'  RL1(1000*( 1-ALPHA)) 

*  4X, 'TRUE  CONFIDENCE  LEVEL* ,5X, ' JACKNIFE' , 

*  /// ,5X, 1 R  SX  N1 ,3X, 'ALPHA  =  . 2' , 12X, ' . 1 ' , 11X, ' . 2' , 

*  12X, 1 .  l' ,/ ,3X,13( ' -' ) ,6X,25( ' -' ) ,4X,21( *  -  * ) , 2X, 16( *  ~  * ) > / ) 

650  FORMAT(3(/),5X,'  CASE  ' ,6X, 1  RL1( 1000*( 1-ALPHA) ) 

*  4X,  TRUE  CONFIDENCE  LEVEL'  ///,5X,’R  SX  N',3X. 

*  'ALPHA  =  .  2'  , 12X . ' .  1 1 , 11X,  .  2 1 , 12X, 1 .  l' ,/ ,3X, 13( ' - ' ) , 6X, 

*  25 (,-'),4X>21(1 -'),/) 

700  F0RMAT(3X,F4.  3,2X,F3.  1,2XJI2,4(7X,F7.  4)) 

710  F0RMAT(3X,F4.  3, 1X,F4. 1,2X, 12 ,4(7X,F7. 4) , 1X,2(2X,F6.  4) ,2(2X,F9.  5)) 
800  F0RMAT(/) 

RETURN 

END 


*5'f.Wf!V*:'f*Vj***********Vf*5V*Vf****Vf***********i'f********Vf****!Wf*******yf*Vf**** 

* 


SUBROUTINE  ZTABL1(ZALPHA,  ALPHA) 

«  SUBROUTINE  ZTABL1  COMPUTE  RIGHT  CUMULATIVE  PROBABILITY  FROM  » 
«  RIGHT  PERCENT  POINT.  » 

REAL  PI,  ALPHA,  ZALPHA 

PARAMETER  (  PI  =  3. 141592  ) 


IF  (  ZALPHA  .EQ.  0.0  )  THEN 
ALPHA  =  0.5 

ELSE 

ALPHA  =  EXP(  -2  *  ZALPHA  **  2  /  PI  ) 

ALPHA  =  SQRT( 1  -  ALPHA  *  (1  +  2*(PI-3)*ZALPHA**4/(3*PI**2))) 
ALPHA  =  0.  5  *  (  1  -  ALPHA  ) 

IF  (  ZALPHA  .LT.  0.0  )  THEN 
ALPHA  =  1  -  ALPHA 
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END  IF 

ENDIF 

RETURN 

END 


*  *  it 

****Vr************************iV******Vf**********VcVf*********5VA**A**Vf****** 

SUBROUTINE  ZTABL2( ALPHA, ZALPHA) 

C 

C  «  SUBROUTINE  ZTABL2  COMPUTES  RIGHT  PERCENT  POINT  ZALPHA  FROM  » 

C  «  RIGHT  CUMULATIVE  PROBABILITY  ALPHA  » 

C 

REAL  ALPHA,  ZALPHA 

IF  ((ALPHA  .GT.  0.0)  .OR.  (ALPHA  .  LT.  1))  THEN 
W  =  -  ALOG( 4  *  ALPHA  *  (  1  -  ALPHA  )) 

ZALPHA  =  SQRT(W  *  (2.06118  -  (  5.72622  /  (W  +  11.6406)))) 

IF  (ALPHA  .GT.  0.5)  THEN 
ZALPHA  =  -  ZALPHA 

ENDIF 

ENDIF 

RETURN 

END 


******5V*5V-.,f*****Vf*:Wc********************************yf******************** 
*  it 

iiititiiitititm,  ■'r?V**********************Vc**Vf****************************A****** 

SUBROUTINE  TTABLE( ALPHA,  NU,  TALPHA) 

C 

C  «  SUBROUTINE  TTABLE  COMPUTES  RIGHT  CUMULATIVE  PROBABILITY  » 

C  «  TALPHA  FROM  T-DISTRIBUTION  BY  GIVEN  ALPHA  AND  NU  » 

C 

REAL  ALPHA,  ZALPHA,  TALPHA,  Al,  A2,  A3 
INTEGER  NU 

IF  (.(ALPHA  .GT.  0.0)  .OR.  (ALPHA  .  LT.  1))  THEN 

CALL  ZTABL2( ALPHA,  ZALPHA) 

Al  =  (ZALPHA**2  +  1)  /  4 

A2  =  (5*ZALPHA**4  +  16*ZALPHA**2  +3)  /  96 

A3  =  (3*ZALPHA**6  +  19*ZALPHA**4  +  17*ZALPHA**2  -  15)  /  384 


TALPHA  =  ZALPHA  *  (1  +  Al/NU  +  A2/NU**2  +  A3/NU**3) 

IF  (  ALPHA  .GT.  0.5  )  THEN 
TALPHA  =  -  TALPHA 

ENDIF 
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ENDIF 


RETURN 

END 


************************************************************************ 
*  * 
************************************************************************ 


SUBROUTINE  SHSORT(A,KEY,N) 
DIMENSION  A(N) ,KEY(N) 

Ml=l 

6  M1=M1*2 

IF  (Ml  .LE.  N)  GO  TO  6 

Ml=Ml/2-l 

MM=MAX0( Ml/2,1) 

GO  TO  21 

20  MM=MM/2 

IF  (MM  .LE.  0)  GO  TO  100 

21  K=N-MM 

22  DO  1  J=1,K 
II=J 

11  IM=I I+MM 

IF  (A(IM)  .GE.  A(II))  GO  TO  1 
TEMP=A( II) 

IT=KEY( II) 

A(II)=A(IM) 

KEY(II)=KEY(IM) 

A(IM)=TEMP 

KEY(IM)=IT 

II=II-MM 

IF  (II  .GT.  0)  GO  TO  11 
1  CONTINUE 
GO  TO  20 
100  RETURN 
END 


*  * 


SUBROUTINE  VAR(Z,  NUMZ,  ZBAR,  ZVAR) 

«  SUBROUTINE  VAR  WILL  COMPUTE  THE  MEAN  AND  THE  VARIANCE  » 

«  OF  ARRAY  Z  >> 

REAL  Z(NUMZ) ,  SUMZ,  SUMZ2,  ZBAR,  ZVAR 

SUMZ  =  0 
SUMZ2  =  0 

DO  100  I  =  1,  NUMZ 

SUMZ  =  SUMZ  +  Z(I) 
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SUMZ2  =  SUMZ2  +  Z(I)  **  2 
100  CONTINUE 

ZBAR  =  SUMZ  /  NUMZ 

ZVAR  =  ABS(  SUMZ2  /  NUMZ  -  ZBAR  **  2  ) 

RETURN 

END 


A-.V****;';***************************************************************** 
*  * 
***5V****yf************A******************Vf*Vc*******AA************?V***yc*iWf 


SUBROUTINE  FINDJ(A,  R,  J) 

«  SUBROUTINE  FINDJ  FINDS  THE  INDEX  OF  ARRAY  A  WHICH  THE  » 

«  VALUE  OF  IT  IS  CLOSEST  TO  R.  » 

REAL  A(1000),  R,  VALUE 
INTEGER  J 

VALUE  =  ABS(A( 1000)  -  R) 

DO  100  I  =  999,  1,  -1 

IF  (ABS(A(I)  -  R)  .LT.  VALUE)  THEN 
VALUE  =  ABS( A( I)  -  R) 

ELSE 

J  =  I  +  1 
RETURN 

END  IF 

100  CONTINUE 

RETURN 

END 


fryf**&*vr*vt****yc*yfyr**yrycVf  ****!&**  yfyfyr**yf  ************************************ 

* 

;V*rt*V?*y?***Vc*******A*A****rtyf*****?V******A*rt*y:ycVcyc****:fc******yf**&**rt*****yr 


SUBROUTINE  TNYVAL( R , S I GMAX , NUMX , ALPHA , Z ALPHA , MUX , YO , 

TBRL1,TCIR1) 

«  SUBROUTINE  TNYVAL  IS  THE  CASE  WHEN  X  IS  TRUNCATED  NORMAL  » 
«  AND  Y  IS  GIVEN  A  VALUE  » 


REAL  R,  SIGMAX,  ALPHA,  TALPHA,  MUX,  P(2),  ZP(2),  A(2) 

REAL  X(100),  XI,  TEMPI,  TEMP2 ,  YO 

REAL  TARLl(lOOO),  KEY(IOOO),  TRLl(lOOO),  TBRL1(2) ,  TCIR1(2) 
INTEGER  NUMX,  CASE,  COUNT 
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DATA  P  /  .95,  .90  / 

DATA  ISEED,  JSEED  /  4875,  7981  / 

DO  400  I  =  1,  2 

PCI)  =  1  -  PCI) 

COUNT  =  0 

CALL  ZTABL2(P(I) ,  ZP(I)) 

A(I)  =  MUX  +  ZP(I)  *  SIGMAX 
DO  200  J  =  1,  1000 

100  CALL  LNORMC ISEED,  XI,  1,  2,  0) 

XTEMP  =  MUX  +  SIGMAX  *  XI 

IF  (XTEMP  .LE.  A(I))  THEN 
COUNT  =  COUNT  +  1 
X( COUNT)  =  XTEMP 

ELSE 

GOTO  100 
ENDIF 

IF  (COUNT  .LT.  NUMX)  GOTO  100 
CALL  VAR(X,  NUMX,  XBAR,  XVAR) 


C  /*  NOW  WE  COMPUTE  THE  LOWER  CONFIDENCE  BOUND 

TEMPI  =  (XBAR  -  YO)  /  SQRT(XVAR*(NUMX-1)  /  NUMX) 
TEMP2  =  1.  /  NUMX  + 

*  ((XBAR  -  YO)  **  2)  /  (2*(NUMX+1)  *  XVAR) 

TEMPI  =  TEMPI  -  Z ALPHA  *  SQRT(TEMP2) 


C  /*  SUBROUTINE  ZTABL1  COMPUTE  RIGHT  CUMULATIVE  PROBA- 

C  /*  BILITY  FROM  RIGHT  PERCENT  POINT. 

CALL  ZTABL1( TEMPI,  TARLl(J)) 

TARLl(J)  =  1.  -  TARLl(J) 

200  CONTINUE 


/*  NON-IMSL  LIBRARY  'SHSORT'  WILL  SORT  TARL1  BY  SHELL 
/*  SORT  ALGORITHM. 

CALL  SHS0RT(TARL1,  KEY,  1000) 

DO  300  J  =  1,  1000 
TRL1(J)  =  TARLl(J) 


noon  *  *  *  oo 


300 


CONTINUE 


C  /*  SUBROUTINE  FINDJ  FINDS  THE  INDEX  OF  TARL1  WHICH  THE  */ 

C  /*  VALUE  OF  IT  IS  CLOSEST  TO  R.  */ 

CALL  FINDJ(TARL1,  R,  J) 


/*  DIVIDING  THE  INDEX  BY.  1000  WILL  GIVE  US  TRUE  */ 
/*  CONFIDENCE  LEVEL.  -  */ 
TCIRl(I)  =  J  /  1000. 


MM  =  1000  *  (1  -  ALPHA) 
TBRLl(I)  =  TRLl(MM) 

400  CONTINUE 

RETURN 

END 


itisititititititititititlsitMcieiritMeMtMekMeMtieMtMcit'MtMtMeMtMeitiekitit  ******************** 

it 

i:itititi;ici<itititi<iti:ititititititititM;itititititt'ti:ititiciticititi<iti;iticititititititicicitititi:ici<Mtitititiritici:icititititic 


SUBROUTINE  JKNIFE(X,  N,  YO,  ZHAT) 

«  THIS  ROUTINE  IS  FOR  THE  'JACKNIFE'  METHOD.  » 

«  DUMMY  PARAMETER  IS  ZHAT.  » 

REAL  X(N) ,X0MIT( 100) ,XOMBAR( 100) ,X0MVAR( 100) ,ZOMHAT( 100) , ZHAT 
INTEGER  M,  N 

M  =  N  -  1 
SUM  =  0. 

DO  200  I  =  1,  N 

DO  100  J  =  1,  M 

IF  (J  .GE.  I)  THEN 

XOMIT(J)  =  X(J+1) 

ELSE 

XOMIT(J)  =  X(J) 

END  IF 

100  CONTINUE 

CALL  VAR(XOMIT,  M,  XOMBAR(I),  XOMVAR(I)) 

ZOMHAT(I)  =  (  XOMBAR(I)  -  YO  )  /  SQRT(XOMVAR(I)) 

SUM  =  SUM  +ZOMHAT(I) 

200  CONTINUE 

ZHAT  =  SUM  /  N 
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RETURN 

END 
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